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Abstract. We consider the dynamics of spin facilitated models of glasses in the non- 
equilibrium aging regime following a sudden quench from high to low temperatures. 
We briefly review known results obtained for the broad class of kinetically constrained 
models, and then present new results for the behaviour of the one-spin facilitated 
Fredrickson- Andersen and East models in various spatial dimensions. The time 
evolution of one-time quantities, such as the energy density, and the detailed properties 
of two-time correlation and response functions are studied using a combination of 
theoretical approaches, including exact mappings of master operators and reductions 
to integrable quantum spin chains, field theory and renormalization group, and 
independent interval and timescale separation methods. The resulting analytical 
predictions are confirmed by means of detailed numerical simulations. The models 
we consider are characterized by trivial static properties, with no finite temperature 
singularities, but they nevertheless display a surprising variety of dynamic behaviour 
during aging, which can be directly related to the existence and growth in time of 
dynamic lengthscales. Well-behaved fluctuation-dissipation ratios can be defined for 
these models, and we study their properties in detail. We confirm in particular the 
existence of negative fluctuation-dissipation ratios for a large number of observables. 
Our results suggest that well-defined violations of fluctuation-dissipation relations, 
of a purely dynamic origin and unrelated to the thermodynamic concept of effective 
temperatures, could in general be present in non-equilibrium glassy materials. 
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1. Why study the aging regime of spin facihtated models? 
1.1. A brief survey of kinetically constrained models 

This paper is concerned with the dynamics of spin facilitated models of glasses in the 
non-equilibrium aging regime following a sudden quench from high temperature to the 
very low temperature glassy regime. Spin facilitated models belong to the broader family 
of kinetically constrained models (KCMs). These are simple statistical mechanics models 
which display many of the dynamical features observed in real glassy materials, such 
as supercooled liquids [Ij, spin glasses [2], or soft disordered materials [3]. KCMs are 
generically defined from a simple, usually non-interacting, Hamiltonian. The complexity 
of glasses is encoded in specific local dynamical rules, or kinetic constraints. For an 
extensive review of early results on KCMs see [1]. 

In this paper we focus on spin facilitated models, in particular Fredrickson- Andersen 
(FA) [5] and East models [6], but we expect that similar behaviour to the one we describe 
here will also be found in other KCMs such as constrained lattice gases. The main 
insight of FA [5] was to devise models that are simplistic, as compared to realistic 
interacting molecular systems, but whose macroscopic behaviour was in agreement 
with the phenomenology of liquids approaching the glass transition [5l [7], displaying 
a super-Arrhenius increase of relaxation timescales on decreasing the temperature and 
non-exponential relaxation functions at equilibrium. Early studies also demonstrated 
that when suddenly quenched to very low temperatures the subsequent non-equilibrium 
aging dynamics of the models compares well with experimental observations on aging 
liquids |H1 [9] . Initially it was suggested that FA models would display finite temperature 
dynamic transitions similar to the one predicted by the mode-coupling theory of 
supercooled liquids but it was soon reahzed that most KCMs do not display such 
singularity, and timescales in fact only diverge in the limit of zero temperature [H [lOl [TT] . 
This implies in particular that after a quench to any non- zero temperature these systems 
are eventually able to reach thermal equilibrium. The time window of the aging regime, 
however, becomes large at low temperatures, and the detailed study of this far from 
equilibrium aging time regime will be the main subject of this paper. 

There has been much interest in KCMs recently. This is partly due to the 
realization [121 [ISl Ull [HI [IS [T71 |18] that their dynamics is characterized by 
dynamic heterogeneity, that is, non-trivial spatio-temporal fiuctuations of the local 
relaxation, which is also a hallmark of supercooled liquids [19]. Many of the studies 
relating to dynamic heterogeneity in KCMs are very recent indeed, having appeared 
since the review [1] was compiled. These studies characterize in great detail the 
heterogeneous dynamics of KCMs. They include: papers defining and quantifying 
relevant dynamic lengthscales in KCMs in equilibrium and their relation to relaxation 
timescales [ini[IIl[2ni[ai[22l[23liai[25l[26l[271[28l^ studies of more quahtative or 
phenomenological consequences of kinetic constraints, including dynamic heterogeneity 
and activated dynamics, for the physics of glassy systems [311 1321 EHl [3ll [35l [36] ; and 
the definition and analysis of new KCMs where kinetic rules, or lattice geometry, are 
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tuned to explore in more detail the range of possible behaviours that can be observed 
in glass models [371 ES [39], UHl [HI [12]. These numerous recent studies have in turn 
been instrumental in encouraging numerical and experimental efforts to measure and 
characterize in more detail dynamic heterogeneity in supercooled liquids [ISj [HI [l5l [l6l 
[13 [13 [13 [50], granular materials [SH [52], colloidal systems [50l [53l [Ml [55] , and soft 
glassy materials [56l [STJ [58], [59] . 

1.2. Aging dynamics in glassy systems: early studies and open questions 

The previous section suggests that dynamic heterogeneity and activated dynamics, 
which have been well-studied and characterized at thermal equilibrium in various 
models and systems approaching the glass transition, play central roles in glassy 
dynamics. However, when moving deeper into the glass phase, glassy materials cannot be 
equilibrated anymore on experimental or numerical timescales. In this non-equilibrium 
state, physical properties are not stationary, and the system displays aging behaviour [2]. 
Experimentally, aging has been well studied at the macroscopic level, in systems as 
diverse as polymers [60], structural and spin glasses [2], and soft materials [3]. A 
full understanding of the non-equilibrium glassy state remains a central theoretical 
challenge [2]. 

Theoretical studies of mean-field models have provided important insights into the 
aging dynamics of both structural and spin glasses [2], [61], [62]. In mean-field models, 
thermal equilibrium is never reached, and aging proceeds by downhill motion in an 
increasingly fiat free energy landscape [63]. Time translational invariance is broken, 
and two-time correlation and response functions depend on both their arguments. 
The fluctuation-dissipation theorem (FDT), which relates equilibrium correlation and 
response functions, does not apply in the aging regime, but a generalized form is shown 
to hold [M] • This is defined in terms of the two-time connected correlation function for 
some generic observable A{t), 

C{t,t^) = {A{t)A{t^)) - {A{t)){A{t^)), (1) 

with t >t^, and the corresponding two-time (impulse) response function 

MAjt)) 
6h{t^) 

Here h denotes the thermodynamically conjugate field to the observable A so that the 
perturbation to the Hamiltonian (or energy function) is 6E = —hA, and angled brackets 
indicate an average over initial conditions and any stochasticity in the dynamics. Note 
that we have absorbed the temperature T in the definition of the response. The 
associated generalized FDT is then 

R{t, t^) = X{t, t^)^C{t, t,), (3) 

with X(t, t„) the so-called fluctuation-dissipation ratio (FDR). At equilibrium, 
correlation and response functions are time translation invariant, depending only on 



R{t,t^)=T- 



(2) 

h=0 
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T = t — t^, and equilibrium FDT imposes that X{t,t„) = 1 at all times. A parametric 
fluctuation-dissipation (FD) plot of the step response or susceptibility 



is then a straight line with unit slope. These simplifications do not occur in non- 
equilibrium systems. But the definition of an FDR through Eq. ([3]) becomes significant 
for aging systems [61],|62]. In mean-field spin glass models the dependence of the FDR 
on both time arguments is only through the correlation function X{t, t^) ~ X{C{t, t^)) 
at large times. This led to the idea that aging systems might be characterized by 
an effective temperature [6l], defined in terms of the FDR, Tcs = T/X. Physically, 
relaxation in glassy systems occurs in well-separated time sectors |62j; it is then easy 
to imagine that each sector could be associated with an effective temperature [65] . 
A thermodynamic interpretation of effective temperatures has also been put forward, 
relating them to the concept of replica symmetry breaking |66j . 

Taken together, these results make the mean-field description of aging very 
appealing, and have set the agenda for a large body of numerical and experimental 
work, as reviewed in [67] . The broader applicability of the mean-field scenario of aging 
dynamics remains unclear, however. While some experiments and simulations seem to 
support the existence of well-behaved effective temperatures [68l [691 ED] , other studies 
also reveal the limits of the mean-field scenario. Experiments have for instance reported 
anomalously large FDT violations associated with intermittent dynamics [TT l [72 l [731174] . 
while theoretical studies of model systems have also found non-monotonic or even 
negative response functions [75l [761 [ZH [ZHl [79] , and ill-defined or observable-dependent 
FDRs [HQ]. In principle, these discrepancies with mean- field predictions are to be 
expected, since there are many systems of physical interest in which the dynamics are 
not of mean-field type, displaying both activated processes and spatial heterogeneity. 
It is thus an important task to understand from the theoretical point of view when the 
mean-field concept of an FDR-related effective temperature remains viable. 

1.3. Aging and dynamic heterogeneity 

Studying theoretically the interplay between relevant dynamic lengthscales and 
thermally activated dynamics in the non-equilibrium regime of disordered materials 
is clearly a challenging task. This problem has been approached in several different 
ways, as we briefly summarize in this subsection. 

Spin glasses represent a particular class of materials which has been studied 
in detail by means of experiment, simulation and theory [2]. For these systems 
mean-field theory provides a set of detailed predictions, in particular for the aging 
properties [62l [81], without accounting specifically for growing dynamic lengthscales 
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or thermally activated processes. More phenomenological approaches, on the other 
hand, directly focus on real space excitations. The resulting predictions [82|, [831 El] 
differ from those of mean-field, and aging in particular comes to be associated with the 
logarithmic growth with time of a dynamic correlation lengthscale. This lengthscale 
is believed to play a crucial role in memory and rejuvenation processes observed 
experimentally [85l [86| \87\ [88] , a link that has been confirmed numerically [891 190| - 
There exist theoretical attempts to extend the mean-field framework to include spatial 
fluctuations and dynamic lengthscales [9T11921 [931 [Ml [951 [96] . 

The situation is somewhat similar for structural glasses. While mean-field theories 
of aging are by now well-established [61j, several alternative perspectives insist that 
dynamic lengthscales and fluctuations should play an important role [HI [311 [971 198] . 
This includes the aging regime. Just as in the spin glass case, extensions of the mean-field 
framework to finite dimensions are feasible in principle [291 [M [991 [TM[TnTl [1021 [T03l[Tni 
but this remains a very difficult task. 

A natural approach is to consider directly systems with local, finite ranged 
interactions. The study of the aging dynamics of KCMs is therefore very promising, 
as these models combine the necessary features of being defined in terms of (effective) 
microscopic degrees of freedom, having local dynamical rules, and displaying thermally 
activated and heterogeneous dynamics. 

Early studies on the aging dynamics of KCMs were mainly dedicated to exploring 
the generic behaviour of FA models after a low temperature quench, or using more 
complex thermal protocols [HI [9] . The idea was to establish KCMs as reasonable glass 
models. The Kob- Andersen lattice gas [105j was then studied in some detail, with results 
that appeared to be in reasonable agreement with mean- field ideas [1061 1107] . The same 
model was also used as a kinetic model to study granular compaction [108111091111011111] . 
Both FA and East models were studied numerically in [11211113] , with somewhat unclear 
results. While two-time correlation functions were shown to exhibit standard scaling 
properties, non-monotonic response functions were found, leading to apparently ill- 
defined FDRs. However, these studies considered non-connected correlation functions, 
making the resulting FDRs of dubious relevance. The one-spin one-dimensional FA 
model, which we study in detail below, was reconsidered in [1141 1115] with a different 
conclusion, this time indicating that FDT is satisfied for this model at all times. The 
East model was also reconsidered in [1161 1117] . with results similar to that found in a 
family of spin plaquette models [118j . which share many similarities with KCMs [26j . 
Non-monotonic response functions were again found at low temperature, which seemed 
to prevent the existence of genuine FDRs. 

It was later realized that although the non-monotonic response functions were real 
and physical [114j . they did not prevent the possibility of defining FDRs with robust 
scaling properties. The confusion arose from the use of an incorrect definition of the 
FDR, whereby the t^-derivative in Eq. ([3]) was replaced, for numerical convenience, by 
a t-derivative. In the graphical representation in terms of an FD plot this corresponds 
to fixing t^ and letting t vary along the curve, rather than the correct reverse procedure 
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where the later time t is fixed and the earher time t„ is varied. As has been 
emphasized several times [801 IHSj I120^ I121j . this procedure can lead to completely 
different behaviour or at least incorrect numerical values of the FDR. Unfortunately, 
the t„-derivative in its definition makes the numerical measurement of the correct 
FDR a very demanding task; equivalently, it is normally difficult to obtain the step 
response xit^tw) as a function of the earlier time t^ with t fixed. However, novel 
numerical methods for accessing linear response in simulations have recently become 
available |12H 11221 1123[ I124j and have made it possible to systematically study FDRs in 
KCMs |125l I126| I127j . The results are, perhaps surprisingly, much simpler to interpret 
physically, as we shall show throughout this paper. 

1.4- Summary of main results and plan of paper 

The models considered in this paper are characterized by trivial static properties with 
no finite temperature singularity, yet they display a surprising variety of dynamic 
behaviour during the aging which can be directly related to the existence and growth 
with time of a dynamic lengthscale. We find in particular that, as in mean-field models, 
well-behaved fiuctuation-dissipation ratios can often be defined for these models from 
violations of the FDT. A major difference between mean-field models and KCMs is 
that FDT violations are only transient in KCMs, and a crossover towards equilibrium 
is always expected. Therefore, there can be no asymptotic connection between FDT 
violations and thermodynamic properties. Similarly, it is not obvious how to connect 
FDRs to effective temperatures in the case of KCMs. In particular the presence of 
negative fluctuation-dissipation ratios for a large number of observables indicates that 
the interpretation of FDT violations in terms of some generalized thermodynamics may 
not be possible. In fact, the aging behaviour of KCMs found in this and related 
work |125l 11261 1127] suggests that well-deflned violations of fluctuation-dissipation, 
of purely dynamic origin and unrelated to the thermodynamic concept of effective 
temperatures, might generically be present in non-equilibrium glassy materials. 

The core material of the paper is divided into the following Sees. [2] and [3] which deal 
with the FA and East models respectively. For both models we flrst give a qualitative 
description of the physical processes governing the aging dynamics. For the FA model we 
then recall (Sec. 12. ip arguments based on an exact mapping to a diffusion-annihilation 
process which demonstrate that the upper critical dimension is (ic = 2 (rather than 
dc = 4). Below this dimension, i.e. in = 1, it turns out that exact results can be 
obtained in an appropriate long-time scaling regime. These are described in Sec. 12. 2[ 
first in general terms and then applied to the local, global and Fourier mode correlation 
and response functions. The latter case (Sec. I2.2.3P includes the other two as opposite 
limiting behaviours, and simulations confirm the predictions across the entire range. In 
Sec. 12. 31 we turn to dimensions d > d^ = 2, where field-theoretic methods can be used to 
calculate the Fourier mode correlation and response functions. Also these predictions 
compare well with simulation results. We conclude the section on the FA model with 
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the derivation of a useful identity for global response functions (Sec. 12.41) that makes the 
latter much more amenable to numerical simulation. The identity requires only that the 
dynamics be given by local spin flips that are modified by a kinetic constraint and so 
should be of general use for numerical studies of KCMs with non-conserved (i.e. Glauber 
spin-flip rather than Kawasaki spin-exchange) dynamics. 

In the analysis of the East model, different techniques come to the fore, in particular 
master equations based on the independent interval nature of the dynamics. These can 
be specifically adapted to find correlation and response for local (Sec. 13.11) and global 
(Sec. 13. 2p observables. Fourier modes turn out to be less useful for the East model; 
instead we consider non-local (distance dependent) quantities (Sec. 13.31) which reveal 
more clearly the effects of the directedness of the constraint, in particular the very 
intricate spatial structure of two-time correlations and responses. 

We conclude in Sec. H] with a summary and discussion of our results. Some 
mathematical background for the FA Fourier mode analysis in d = 1, and technical 
details of the more involved East model calculations, can be found in the appendices. 

2. Fredrickson- Andersen model 

The FA model [5l [7] describes the dynamics of binary spin-variables nj = 0, 1 on 
a hypercubic lattice in d dimensions. An up-spin (n^ = 1) represents a highly mobile 
region or, because of the associated energy cost, a "defect". Down-spins (nj = 0) model 
dense, immobile regions and are energetically preferred. The energy function is assumed 
to be non- interacting, E = ^ • n^. In a continuous-time Markov dynamics, Glauber rates 
for flipping spin i are prescribed as 



where n = (ni, . . . ,nAr) and c = 1/(1 + e^) is the equilibrium up-spin concentration; 
(3 = 1/T denotes the inverse temperature as usual. Without the factor fi{n) the model 
would consist of non-interacting spins, with rates c and 1 — c for up- and down-flips, 
respectively. All the interesting physics is therefore in the facilitation factor fi{n). It 
is chosen to be independent of Ui and this ensures that detailed balance with respect 
to the energy function E is maintained. Physically, fi represents the assumption that a 
change of state in region i is possible only if some neighbouring regions are in a mobile 
state. We will focus exclusively on one-spin facilitated FA models below, which are 
rather simpler to understand than those imposing facilitation by more than one spin [3]. 
For simulation purposes it is then simplest to set /j(n) = if = for all nearest 
neighbours j of site and otherwise fi{n) = 1 [128] . For analytical work it is more 
convenient to use the form of fi suggested originally O [7] , 



j=n.n.(i) 

which is again zero when all n.n. sites are in immobile states, but has nonzero rates 
increasing in proportion to the number of mobile neighbours. Physically, the key 



= /i(n)[(l - c)ni + c(l - n^)]. 



(6) 




(7) 
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influence of fi is to make moves witliout mobile neighbours impossible; the difference 
between the two above definitions of fi for the rates of the allowed moves therefore has 
no qualitative effects on the observed behaviour. 

From the point of view of glass modelling the interesting parameter regime of FA 
models is that of low c, where mobile regions are few and far between. After a quench 
from high temperature (corresponding to initial defect concentration cq = c(T oo) = 
1/2) any such defects that are not isolated from each other will quickly flip down, and 
a state is reached where essentially all defects are isolated from each other. From then 
on, any further reduction in the defect concentration 



proceeds via a diffusion-coagulation process. Diffusion of defects can take place when 
a defect is created next to an existing one, with rate c. Both defects are now mobile 
and can flip down quickly (rate 1 — c ~ 1); if the original defect does so first, it leaves 
the other one behind and has effectively moved by one site. The diffusion constant for 
this process is D = c/2, with the factor 1/2 accounting for the probability of the newly 
created defect flipping back down first. When two defects diffusing by this mechanism 
meet, one of them can flip down, leading to defect coagulation. The reverse process of 
branching is also possible, by creating two defects next to an existing one which then 
flips down. However, as two new defects are involved the effective rate is ~ and 
therefore negligible compared to those for diffusion and coagulation at least in the time 
window c^^ ^ t ^ c~^. For larger times branching must of course eventually become 
important since it is the only process that creates defects and so is able to stabilize the 
defect density at its equilibrium value, n(t — > oo) = c. 

Conceptually, it is important to note that the effective diffusion-coagulation 
dynamics has a dynamic critical point at c ^ 0, where both timescales and lengthscales 
diverge [211 [23]. One can define appropriate critical exponents; this is simplest from the 
equilibrium dynamics at low c. The order parameter exponent n(t oo) = c ~ 
is fixed to (3 = 1 because of detailed balance. The lengthscale ^ and timescale 
trei of appropriate relaxation functions define the other exponents as ^ ~ c"'^ and 
Dtj-ci ~ ~ c~^'^. We note as an aside that, in the context of dynamic criticality, 
the asymptotic FDR X°°, which is defined as the limit of X{t,t^) for widely separated 
times t 3> tw, appears as a universal amplitude ratio that can distinguish different 
dynamic universality classes [129] . 

2.1. Mapping to diffusion- annihilation 

Initial field-theoretical studies [23] suggested that the FA model was in the dynamic 
universality class of directed percolation, with an upper critical dimension dc = 4. 
Closer inspection reveals, however, that there is a hidden symmetry which reduces this 
value to dc = 2 [23]. This can be deduced from an exact mapping between the FA 
model and another defect-diffusion model where defects or "particles" annihilate in 




(8) 
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pairs {A + A 0, where A symbolizes a defect) and are similarly created {0 ^ A + A) 
in pairs. 

The mapping is most easily derived in a quantum mechanical representation of 
the Markov dynamics |130l 11311 1132] . One associates to each possible state n of the 
system a distinct unit vector \n) in a vector space, and to the time- dependent probability 
distribution over states pt(n) the vector \pt) = J2nPt{''^)\'^) ■ (The basis vectors |n) have 
unit length and are mutually orthogonal.) It is then easy to check that the dynamical 
evolution of this vector takes the form dt\pt) = W\pt) with the master operator 

W = - - C)n. + C(l - (9) 

i 

Here, = X]ti^*I^)(^I operator measuring the value of the local spin, 

fi = Sj=n n (i) "-i' Operator that flips spin i, Fi\ni, . . . , rij, . . . , wat) = 

\ni, . . . ,1 — Hi, . . . , n^v). Averages are obtained by multiplying with the uniform state 
{^\ = Z^ri("'l fro"^ the left, e.g. {ni{t)) = {e\ni\pt). 

The key mathematical statement that establishes the mapping mentioned above is 
that there is an operator V such that the similarity transformation VWV~^ of the FA 
master operator produces the master operator of an annihilation-pair creation process 
with appropriate rates [2lj. One can then verify directly that the equilibrium correlation 
functions (more precisely those involving only two separate times) of the two models 
are identical up to trivial factors. As a consequence, also the dynamic exponents and 
dynamic universality class of the two models must be the same. In particular, one has 
dc = 2 and the critical exponents are [z, v,l3) = (2, l/c?, 1) in < 2 while they take the 
mean-field values {z, Vil3) = (2, 1/2, 1) in higher dimensions. 

The feature of the annihilation-pair creation process that is responsible for the 
upper critical dimension being two rather than four is the conservation of parity: if the 
number of particles is even to start with it remains so for all times, and similarly if it is 
initially odd. Mathematically, the parity operator, 

P=(-l)E.n.^ (10) 

commutes with the master operator, and maps two different steady states onto each 
other (which differ only in the sign of the probability amplitude of states with an 
odd number of particles). The similarity relation between annihilation-pair creation 
and the FA model implies that the latter also has a symmetry; this is the "hidden" 
symmetry mentioned above. Its physical nature is somewhat difficult to elucidate since 
the symmetry operator that commutes with the FA master operator is not diagonal 
and so does not correspond to a standard observable obeying a conservation law. 
Mathematically, the symmetry swaps terms in the master operator across n.n. bonds: 
the term describing a flip of spin i facilitated from a neighbouring site j is mapped to 
the one encoding flips of spin j facilitated from site i. This insight makes clear that 
the symmetry also applies to the East model discussed below, except that the mapping 
here converts the East model into the West model, where the direction of facilitation is 
reversed. While this does not lead to obvious ways of calculating the correlation and 
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response functions that we will study below, it does provide some general restrictions. 
For example, general multi-point two-time equilibrium correlations 

r s 

(nKw-c]n[^>(o)-^])' (11) 

p=l (T=l 

can be shown to vanish unless the rightmost spins appearing at time and t are identical, 
i.e. unless maxp(ip) = ma.x„{j„). For two-point correlations (r = s = 1) this reduces to 
the well-known fact that only local correlations are non-vanishing at equilibrium. 

The above mapping can be extended in a number of ways, for example to 
generic coagulation-branching and annihilation-pair creation models, and to "bosonic" 
processes where each site can be occupied by more than one particle |2l]. For the 
original "fermionic" (hardcore repulsion) versions the mapping also has a very appealing 
geometric structure, with all operations being effectively rotations on a "spin sphere". 



Rr{t, tw) = T 



(13) 

hi=0 



2.2. Exact results in d = 1 

In this section we present exact long-time scaling forms for FDT violations in the Id 
FA model. More precisely, the dynamics following a quench from a random initial state 
to low temperature c <C 1 is considered, in the nonequilibrium aging regime of times 
c^^ ^ t ^ c^^. To probe violations of FDT we introduce the connected correlation and 
response functions, 

Cr(t,t^) = {ni+r{t)ni{t^)) -n{t)n{t^), (12) 
5hi{t^) 

where temperature is again included in the definition of the response function. The fields 
hi are conjugate to the n,, so that in their presence the energy becomes E = nj(l— /ij). 
We will also consider the step-response Xr(t,t^) = dt' Rr(t,t') to a perturbation 
applied from time tw onwards. 

Our analysis is based on the fact that defects effectively diffuse and coagulate on 
the 0{c~^) timescale. As explained above, the rate for diffusion is -D = c/2 while 
coagulations occur with rate 7 ~ 1. We will, from now on and throughout this section, 
measure time in units of (except when stating simulation parameters) so that -D = | 
and 7 ~ ^ 1. The coarsening dynamics of this diffusion-coagulation process at long 
times t 1 becomes independent of 7 |133[ 1134] : defects typically first have to diffuse 
for a time 0{t) before they occupy adjacent sites where coagulation is possible. The 
duration of this reaction, which depends on the rate 7, is only 0{1) and therefore 
negligible in comparison: the process is said to be in the diffusion controlled regime 
|133j . The irrelevance of the precise value of 7 at long times allows us to adjust its value 
to 7 = 1, the reason for this choice being that (only) the process with D = j admits an 
exact solution |131l I132[ 1135] . Our results will then be exact to leading order at large 
time t ^ 1 for any process with 7 > 0; this includes the FA model. Of course, the 
diffusion-coagulation picture applies to the FA model only so long as t <^ is still 
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satisfied. In the FA dynamics, branching events become important on this times scale 
and initiate a crossover to equihbrium. Therefore, whenever we talk about long time 
behaviour for t — oo it is understood that the low temperature limit c — is taken first 
(in order to remain in the regime 1 -C t -C c~^; recall that our time-unit is now c~^). 
Throughout this section we will use the symbol "~" to indicate results which become 
exact in this long-time limit. 

The scaling of the connected two-time correlation Cr{t,t„) in the effective FA 
process with D = 7 = | was derived in |127j . The result reads 

00 

i,j=0 

00 

- J2 Gi+r,+r-i),ii,j){r) e-'*-[/. + /i+i](2t„) [6j,o + H,{2t^)] 

i,j=0 
00 

-2t„ 



E 



G(_,,_,_i),(i,,-)(r) e-^*-[J, + /i+i](2t„) [5,-o + H,{2t^)] , (14) 

i,j=0 

where In{t) denotes modified Bessel functions, Eq. ( lA.ip . and the explicit form of the 



functions Hn{t) is given in Eq. (]A.2p . To save space we use the short-hand notation 
[■](x) to indicate that all functions contained in the square brackets have the same 
argument x. The Gij{T) with i = (21,^2) and j = (ji, j2) are Green's functions. 

The result Eq. ( fT4l) applies to the FA model to leading order in ^ 1 but uniformly 
in r > 0, i.e. including small values of r. 

An expression for the response function Rr{t,t^) is also given in [127] . Before 
stating the result we consider the effect of the associated perturbation. It is convenient 
to express the response function in the operator formalism from above, Rr{t,t^) = 
{e\ni^re^'^Vi\pt^) , where Vi = TdW/dhi\hi=o and W the master operator of the effective 
FA process. The perturbation operator Vi accounts for the fact that the field hi reduces 
the energy barrier for creation of a defect at site i; recall that E = ^^^^(l — hi) in 
the presence of the perturbation. Consequently the diffusion rates for entering site i are 
increased. The rates for leaving site i, however, are not affected by hi because they relate 
to processes involving the creation of a defect at sites i — 1 or i + 1. This locally directed 
perturbation can be represented as the sum of two mechanisms |127] : (i) an asymmetric 

(a) 

perturbation that increases the rates for entering site i but decreases the ones for 
leaving it and thus traps diffusing defects, and (ii) a symmetric perturbation l^'-'^'' which 
increases both the rates for entering and leaving site i and thus accelerates the local 

(a) (s) 

dynamics. In combination Vi = reproduces the directed perturbation of the 

local diffusion rates via the field hi] the two contributions can be written explicitly as 

^/'^ = \iFi_,F, - l)(n,_i - hi) + ^(F,+iFi - l)(n,+i - n,), (16) 
Vl'^ = i(F,_iF, - l)(n,_i - h,)' + i(F,+iF, - l)(n,+i - h,)^. (17) 
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Clearly an analogous decomposition then applies to the response function 

Rrit,t„) = Ri^\t,t^) + Rl^\t,t^), (18) 

with Ri^^^\t,tvj) = {e\ni+re^'^v}^^^^ \pt^) the contributions from the asymmetric and 
symmetric perturbations. We will see below that Ri'^\t,t^) and Rl^\t,t^) have very 
different scaling properties, each being dominant in a different time regime. For now we 
simply state the relevant expressions [127] . 

Ri^\t, t^) ~ dt„^e-''lr{t - tw)[/r-l + 2/, + /,+i](t + tw), (19) 

Ri'\t,t^) ~ {Ir{t - tw)[-/r-2 + 2Ir " J.+s] (t + ^w) 

+ [Ir-l - Ir+l]it - tw)[/._l - /.+l](t + tw)} . (20) 

Like Eq. (fT^ . Eq. (fT9l) applies to leading order in ^ 1 and uniformly in r > 0. The 
symmetric contribution Ri''\t,ty^), Eq. (1201) . is sensitive to the value of 7 and for this 
reason only applies when both r,t^ ^ 1 [127j . 

2.2.1. Local correlation and response The violation of FDT associated with the local 
defect observable rzj has been considered in all simulation studies of the FA model 
|112[ 11141 I115[ I125j . However, as we now explain, this observable is ill suited to 
the measurement of FDT violations. One can derive the scaling of the relevant 
autocorrelation Co(t,tw) and response Ro{t,t^) functions from Eqs. (HM . ( fTOil and (l20ll 
evaluated at r = 0. In the quasi-equilibrium regime of small r > and large t^ ^ 1 
the autocorrelation scales as 



Co(t,tw) ~ ^^e-Vo(r) ~ n{tMr). (21) 



Here we have identified the defect concentration n{t^) = l/yvrt^ and the random walk 
return probability Pri^) = e~'^/o(r). In the quasi-equilibrium regime defects are typically 
too far from each other to meet and therefore may be treated as independent random 
walkers. Correspondingly, the average {ni{t)ni{t„)) is given by the probability n(t„) of 
having a defect at site i at time t^ multiplied by the probability Prir) for this defect to 
return and occupy the same site at time t. This is then also the leading contribution 
to the local autocorrelation function Co(t,tw); the subtraction of n(t)n(t^) only gives a 
subleading correction because n(t) -C 1. 

In the aging regime of large ^ 1 with fixed r/t^, an expansion of Eq. ffT4l) yields 
the scaling behaviour |127j 

Coit,t^) -^^ I dx I rfy(x + y)e-(^*-/-)(^"+^')$(a; + y) 



JO 

00 POO 



7r3/2 ^2 

r<00 



dx / rf?/(a;-?/)e-(2W-)(x-^+y^)e-^'$(y). (22) 
Jo 



Here ^{z) = (2/^/11) due " denotes the complementary error function. The inte- 
grals in Eq. (!22|) can be solved but the result is not particularly simple; we only state 
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the resulting scaling in the two limiting cases 

r < : Co{t, t„) ^ ^==, (23) 

IT y/ZTt^ 

r»t^: Co(t,U^^^^. (24) 

This completely characterizes the defect autocorrelation Co(t, t„): it has the initial value 
Coitw, tw) ~ '^(^w) from which it decreases according to Eq. fl2T|) . When 1 <^ r <^ t^ this 
quasi-equilibrium result matches with the expansion Eq. ( l23l) . As r is further increased 
beyond t^ the scaling of Co(t,tw) finally crosses over from Eq. (l23l) to Eq. fl2^ . 

The scaling of the defect response function Ro{t,t„) = R^\t,t^) + i?Q^'*(t,tw) is 
unusual in several ways. It is instructive to first discuss the step-response function 
Xo(^)^w)- The latter is always dominated by contributions from the asymmetric 
perturbation Xo(^5^w) ~ X^o\tjtw) as can be shown [127] from Eqs. ( |T9l 1201) . An 
expression for Xo'^''(^; ^^w) = Jt dt' R'^\t,t') is straightforwardly obtained from Eq. (fT9|) . 
In the quasi-equilibrium regime of small r > and large tw ^ 1 one finds the scaling 

xt\t,t^)-^n{t^)[l-pr{r)], (25) 

while Xo^itytw) ~ n{t) in the aging regime ty^ ^ 1 at fixed r/t^- This non-monotonic 
behaviour of Xo^''(t, tw) is due to a simple mechanism: when applying the asymmetric 
perturbation - which traps diffusing defects - at site i and from time t„ onwards, this 
increases the probability {ni{t)) to find a defect at site i. The step response increases 
on a r = 0{1) timescale, c.f. Eq. (l25ll . from zero to n{tv,)- But in the aging regime 
and for r ^ the density of defects n(t) in the system decreases and along with it 
so does Xo\t,t„) ~ n{t). While this makes perfect sense, one has to be careful when 
recovering the response R'^\t,ty^) = —dt^Xo\t^t^) from this scaling result: one would 
incorrectly conclude that R'^\t,t^) = —dt^n{t) = 0. The resolution of the apparent 
paradox is that Xo'^''(^)^w) is dominated by small r contributions of i?o'^^(t, t„), so that 
the aging behaviour of R'^\t,ty^) only shows up in subdominant terms in Xo'^''(^! ^w)- 
Summarizing for the step response, Xoit^tw) ~ Xo'^''('^5 ^^w) is dominated by the trapping 
effect of the asymmetric perturbation and X^\t,ty^), in turn, is dominated by short 
time contributions from R^\t,tv,)- This inconvenient insensitivity to aging effects 
can be avoided by considering directly the response function RQ{t,t^), rather than the 
step response xoit^tw)- From Eqs. (fT9l [20l) its constituent asymmetric and symmetric 
contributions scale as 

^■'(M.O ~ (27) 

in the aging regime T,t^ ^ 1 and r/t^ = 0{1). Note the negative sign of RQ\t,t^). 
To understand this, recall that applying the symmetric perturbation at site i increases 
the local diffusion rates. This in itself does not affect the defect concentration {ni{t)) at 
site i. However, the locally accelerated dynamics increase the likelihood for coagulation 
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with neighbouring defects and thus reduce the local concentration {uiit)). Consequently 
Ro\t,t„) is negative. Considering more quantitatively the asymmetric and symmetric 
response functions, Eqs. ( l26l 1271) . we see that they are of the same order 
when r/t^ = 0{1). For much smaller r <C t„, the trapping effect encoded in 
= C(r-3/2t^-i/2) dominates over i?f^(t,tw) ~ 0{t-^/H^-^/^), and this is 
at the origin of the analogous dominance that we found in the step response. In the 
opposite regime r ^ t„ the situation is reversed, with effects from locally accelerated 
dynamics described by RQ\t,t„) = 0{t^'^) dominating over Rlf^{t,ty^) = C(r^^tw)- 

Based on the above insights regarding the scaling of Co(t,t^) and Ro(t,t„) we are 
now in a position to discuss the resulting violations of the FDT and the usefulness of FD 
plots in revealing them. In the quasi-equilibrium regime of small r > and large t„ 3> 1 
we have that Co{t,t^) ~ n{t^)pr{T), Eq. (I2l]), and Xo{t,t„) ~ xt\t'>t^) ~ '"■(^w)[l - 
Pri^)], Eq. ( l25l) . This satisfies the equilibrium FDT xo(^, ^w) ~ Co{t,t) — Co(t, tw) 
to leading order in t„ ^ 1. Furthermore, the autocorrelation Co(t,t^) decreases 
like Pr(T) = 0(t^^/^) with increasing r and thus drops to a fraction of 
its equal time value Co{t^,t^) by the time we reach the aging regime where r and 
tw are of the same order; the susceptibility Xoit^tw) approaches n(t) up to a small 
deviation of the same order. Consequently FD plots for the local defect observable 
show effectively only the quasi-equilibrium behaviour of Co(t, t^) and Xo(^) ^w), with the 
nontrivial aging regime essentially invisible at long times |114l 11151 I125j . This type 
of scaling makes the local defect observable a poor choice for probing FDT violations, 
especially using FD plots. Regardless of these caveats, however, there is a well defined 
FDR Xo(t,tw) = Roit,t^)/dt^Co{t,t^). It scales as Xq ~ Xo{t/t„) in the aging limit 
[1271 1125j . crossing over from quasi-equilibrium Xq ~ 1 for r <^ t„ to the asymptotic 
FDR Xq ~ X^ for r ^ t^- From the expansions Eq. (|24|) and Eq. (|27|) one obtains the 
negative value 

= ^ -3.307. (28) 

° 67r-16 ^ ' 

We show simulation results in Fig. [1] for various times t. These confirm the findings 

of this section. In particular, the FD plot closely follows the equilibrium FDT except in 

the regime of large time differences, and the non-equilibrium part shrinks into the top 

right corner as t increases. Very detailed measurements in this region, as shown in the 

inset, are consistent with the predicted asymptotic value of the FDR from Eq. ( l28i) . 



2.2.2. Global correlation and response A much clearer picture of the violation of the 
FDT emerges when considering a global observable like the energy [125[ I127[ I136[ 
I137j . We introduce the corresponding normalized connected correlation and response 
functions 



CE(t, t^ 
RE(t, 



- mt)E{t^)) 

T 6{E{t)) 



{E{t)){E{t^))) 



N 5hE{t^ 



(29) 
(30) 
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Figure 1. Normalized FD plots for the local observable in the d ~ 1 FA model at 
temperature T — 0.1 and for times t = 2 x 10^, 5 x 10"*, 10^, and 10^ (from bottom to 
top). The dashed diagonal line indicates the equilibrium FDT. Inset: final part of the 
< = 10® data, enlarged to show the nonmonotonic response and the asymptotic FDR, 
— —3.307 (solid line). From Ref. [125] . Copyright American Physical Society. 



The field He uniformly shifts the energy according to E = (1 — hE)Ylii'''^i- When 
substituting E = Yli Eq. ( !29l) and using translational invariance one readily 

shows that in fact CE(t,t^) = J2r^r(t,t^) and similarly RE(t,tw) = J2r ^r(t,tw) from 
Eq. fl30|) . The scalings of energy correlation and response functions can thus be derived 
by summing Eqs. dH]), (fT9i) and ( l20l) over r. 

It can be shown [127j from Eq. (fT^ that the scaling of energy correlations in the 
aging regime tw ^ 1 with fixed r/t^ is given by 

CE{t,t^)--^l dx I di/(x + 2/)e-(*-/^)(^+^)'$(x + 2/) 



TT T 



, dx / d?/(x-?/)e-(*"/")(^-^)'e-^'$(?/). (31) 

Again these integrals can be solved exactly. The result is bulky but would be required 
in full only if we wanted to inquire into the precise form of the crossover between the 
simple limit behaviours 

r < tw : CE{t, tw) ^ ^ , (32) 

T ^ t„ : CEit,t^) ^ — — -J- — -J-. (33) 

Energy correlations have the plateau value Eq. (!32l) for < r ^ tw- This is because 
defects typically do not meet when r ^ tw and therefore CE(t,t^) ~ CE{ty^,t^). Only 
once r increases to become comparable to tw do coagulation events start to decorrelate 
E{t) from E{tyj), leading to the crossover from Eq. fl32|) to fl33l) . 

The scaling of the energy response function i?£;(t,tw) = R^^\t,ty^) + R^^\t,t^) can 
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also be understood in very simple terms. The response R^^\t,t^) = vanishes exactly 
because, when applied to all sites, the asymmetric perturbations of the diffusion rates 
cancel. Therefore RE{t,t„) = R^^\t,tyj) = Ri^\t,t^). From Eq. (l20l) this gives 

i?s(t,tw) ~ - ^^^3/, ~ dtn{t), (34) 

in the aging regime ^ 1 with fixed r/t^. Due to the normalisation by Rsit, t„) 
measures changes in the defect concentration n{t) = {l/N){E{t)) . Now, in the presence 
of the perturbation all diffusion rates are increased. This is equivalent to giving the 
system some additional time 6t to evolve. The energy response function is therefore 
RE{t,tw) ~ [n'{t + 5t) — n{t)]l5t ~ dtn{t). To understand the sign of this, note that 
the field He decreases the energy barriers for defect creation. While in equilibrium this 
would increase the density of defects, in the nonequilibrium coarsening dynamics this 
speeds up the relaxation of the system - it decreases n{t) and hence REit,t^) < is 
negative throughout. 

From the scalings Eq. (!3T]) and Eq. one obtains the FD plot for the energy. 
Contrary to the case of the local defect observable Ui the energy produces a nontrivial 
FD plot (contained in Fig. [2] below). This can be measured rather easily in simulations 
as discussed in Sec. 12.41 and makes the energy a good observable for probing violations of 
FDT in the FA model. From the scalings in Eqs. (133 1IM1) one shows that the asymptotic 
FDR for the energy is 

XS? = (35) 

It matches exactly the value X^, Eq. fl28l) . obtained from the local defect observable 
but is much more accessible to measurement [12511136] . We will discuss this and related 
points further in the following section. 

2.2.3. Fourier mode correlation and response The above results describe the two 
extreme cases of spatially local and global observables. We will now consider observables 
with a finite intrinsic lengthscale in order to reveal in more detail the physics underlying 
the violation of FDT in the Id FA model. Specifically, we consider the Fourier modes 
of the arrangement of defects, Ug = The associated connected correlation 

and response functions are 

C,it,t„) = ^ [(n,(t)n_,(t„)) - (n,(t))(n_,(t„))] , (36) 



Rq{t, t 



T 5{n,{t)) 



N 6h.g{t^ 



(37) 

The fields hg are conjugate to rig, resulting in a perturbation SE = —hgUg. Clearly, 
Cg(t,t^) is the dynamic structure factor and Rg{t,t„) the conjugate response function. 
Using translational invariance one readily verifies that Cg(t,t^) = J2r^~^'^^^r(^^^^) 
the Fourier transform of the spatial correlations, Eq. ( |T2l) . and likewise for Rq{t, t„). The 
general results Eqs. (fT4|) . (fT9|) and (l20l) are therefore once again the key for analysing 
the functions Cg{t,t^) and Rg(t,t^). 
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In order to derive Cq{t,t^) from Eq. (IT^ the Fourier transforms of the Green's 
functions are needed; these contain all r dependence in Cr{t,t„). Using the identity 
Eq. ( 1A.7I) one immediately has 

r 

G(-)(r;z,j) = ^e-''"-G(.,._i),(,,)(r) = e-'^(^+^+i)^e-2-[J,_,_i - (2rcos|) , 

r 

and in terms of these functions 

oo 

C,(t,t„) ~ e-2*-[/o + /i](2t„) ^ Re [GW(r;z,j)] i^.+,+i(2tw) 

i,j=0 

oo 

- J] 2Re K-nr;^, j)] e-2*"[/, + /,+i](2t„) [5,-o + i^,(2t„)] . (38) 

i,j=0 

In the first line of Eq. fl38l) the imaginary part of G^ci^\T]i,j) drops out since it is 
odd in i — j. The second line, on the other hand, is a sum of two contributions 
{T;i,j) from Eq. ([Til), which is again real. The scaling behaviour of 
Eq. (!38|) in the aging limit t„ — > oo at fixed r/t^ and tw?^ follows from the asymptotic 
expansions Eqs. ( 1A.31 IA.4I) of the functions /„(t) and Hn(t), and is given by 

4 _1t-o2 



rfx / dy{x + y) cos (v^g(a; - y)) e"^*"/^)^''^^) $(x + y) 



dx / c/y (x - ?/) cos (v/t^g(a; + ?/)) e (*»'/^)(^ 2')''e "^^^{y) 



(39) 

Time arguments = Cq{t^t^) were omitted to save space. Energy correlations are 
obtained as the special case CE{t,t^) = Cq=o{t,t^) as is clear from the definition of 
Cq{t,t„), Eq. fl36|) . and indeed Eq. fl39l) with g = reduces to Eq. fl3T|) . Moreover, the 
inverse Fourier transform Cr.=o(t,tw) = /('^Q'/2vr) Cq(t, tw) yields the autocorrelation. 
The resulting scaling of Cr=o(t,t^), Eq. (l22l) . is then also a special case of Eq. (139|1 . 

Equation (!39l) can be rearranged into a more useful form. We rotate the integration 
coordinates u = y + x and v = y — x, whereby the integration domain becomes m G [0, oo] 
and V G [—u, u]. In the first line of Eq. (I39l) the integration over v can then be carried 
out. In the second line, we integrate by parts in u. This leads to some simplifications 
and produces, after combining the v G [—u, 0] and v G [0, u] integration ranges. 



y/nt-^/^ TTT^I^ Jq y/t^q ' 



where (40) 
/(«, .) = !i±^e-K"+^)^$ ( "t^] - !i^e-^(-)^$ ( . (41) 



This is our general scaling result for the dynamic structure factor. Let us now consider 
some limit cases of Eq. (HOj) . We first study the regime r <^ tw In this case the 
exponential e"*^*"^/"^^^^ peaks sharply at f = 0. To leading order in r/t^ we may replace 
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the upper integration limit u of the integration with oo and approximate f{u,v) by 
its expansion at f = 0. Since f{u,0) = we use f{u,v) ~ v dyf{u,v)\y=o. Also noting 
that f{u,v) satisfies the identity dJ{u,v)\^=o = 9„[Me-"'/^$(M/2) - {2/^/^)e~^'/^] and 
using integration by parts in u then gives 



e 2^'''' - A/2e"^*"'^' + J du cos (y/Uqu) u e (^) 



(42) 



in the regime t <^ t^ and to leading order in r/t^- It turns out that Eq. (I42l) reduces 
correctly to the static structure factor for r — > 0; this can be verified rigorously by 
setting r = in the full expression Eq. (1551) and then expanding for large t^ and fixed 
t^q"^. Either way one ends up with the result 

1 - v^e"^*"'''' + j du cos (VU^Qu) ue'^""^^ Q 



Cqit^f,, t 



W5 J 



V Trtw . Jo ^ ' ^ 2 ^ _ 

The static structure factor has a local minimum at g = 0, where it measures energy 
fiuctuations C£;(tw,i^w) = Cq=Q{t^,t^). We find the value Cg=o(^w, ^w) = (3-2V2)/V7rt„ 
from Eq. (H3il . This is, of course, consistent with the expansion Eq. (1321) of CEit,t^) 
from above. For t^q"^ ^ 1, Cq(t^,t^) quickly approaches n(t^) = Ij ^iit^. This large 
q behaviour originates from the scaling of Cr=Q{t^,t^) = n(t^) — nit^Y ~ nit^). The 
local minimum in the static structure factor is caused by an "effective repulsion" of 
defects in the diffusion-coagulation dynamics: defects that have survived up to time 
are likely to be far from each other, i.e. at a distance 0{^/tZ). A similar effect was seen 
in the aging regime of plaquette models |126j . 

To understand the scaling of Cq(t,t^) for 1 ^ r ^ t„ we return to Eq. ( l42l) . The 
profile of Cq(t,t^) is as follows: on the scale t^q"^ = 0{1), where rg^ ^ 1 since r <^ t^, 
Eq. fH2]) is still essentially equal to the static structure factor Eq. fH5]l . However, on the 
larger scale in q where rg^ = 0{1), and hence t^q"^ ^ 1, one has from Eq. fH2l) . 



C,(t,t„)^^e-^«l (44) 

This makes sense intuitively. The wavevector scale rg^ = 0{1) corresponds to a 
length scale £(r) = 0{^/T). Here, since r ^ t„, we may treat defects as independent 
random walkers. Correlations in Cr(t,t„) on the lengthscale £{t) are diffusive, and 
this is consistent with Eq. (l44l) . On the lengthscale £(tw) = 0{-\/t^), corresponding to 
twg^ = on the other hand, spatial correlations Cr{t,tw) are still essentially static 

and this translates into Cq{t,t^) ^ Cq(t^,t^), c.f. Eq. P2l) and Eq. fH3l) . 

We now turn to the scaling of Cq{t,t^) in the opposite regime r ^ t„. Here 
an expansion of Eq. (HOl) follows rather straightforwardly. The exponential e"^*"/"^^^^ 
becomes flat and may be replaced by unity to leading order in t^/r; this is justifled 
since f{u, v) vanishes sufficiently fast in v. Moreover, the only relevant wavevector scale 
is rg^ = C(l) so that sm{y/t^qu) ~ ^/t^qu. Here we have used ^/tZq -C 1 since t^ <^t 
and the fact that again only small u contribute to the integral. Thus altogether. 



Cq{t, tv 



^3/2 



oo 



duu dvvf{u,v) 



(45) 
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The expression in the square brackets only produces an overall coefficient. To evaluate 
it one has to change the integration variables back to x, y. Then the integrals can be 
factorized and one readily obtains 

C,(i,U«|^^e-H', (46) 



in the regime r ^ t^ and to leading order in t^/r. This result contrasts with Eq. (I44ll 
in several ways. First, the scaling of the amplitude of correlations Cg(t,t„) changes 
from (9(tw^^^) in the regime r -C t„ to 0{t^/T^/'^) when r ^ t^. Second, the shape of 
Cq{t, t„) on the wavevector scale rq^ = 0(1) is always Gaussian, however, the exponents 
in Eq. (H4l) and Eq. (H6ll differ by a factor of |. Third, the overall coefficient in Eq. 
is noi related to the static structure factor at g = 0. These differences between Eq. 
and Eq. fH6|) arise from the fact that for r ^ many-body effects play an important 
role; defects meet and coagulate. This regime is fluctuation-dominated with all loop 
diagrams contributing in a field theoretic framework, see Sec. 12.3.41 below. 

Having discussed the dynamic structure factor Cq{t,t^) in detail we now turn to 
the analysis of response functions. Let us ffist consider the asymmetric part Bi^\t,t^) 
of the response function. The Fourier transform of Eq. (fT9l) follows immediately from 
the identity Eq. flA.Sp and is given by 



-2t 



(47) 



where A = ^Jt'^ cos(g/2)^ -|- sin(g/2)2. The scaling of this expression in the aging 
limit t^ ^ (X) with t jt^ and t^q^ fixed follows using Eq. (1A.3P as 

B^it^t^ ~ 5t„^e-i(^-*wA^)V. (48) 



At g = this reduces to the asymmetric part of the energy response R^^\t,t^) = 
-R^^o(^;^w) and vanishes, c.f. Eq. (HHj) . since - as discussed above - the asymmetric 
perturbations cancel when applied uniformly. Integration over q, on the other hand, 
yields the scaling Eq. fl26l) for the auto-response Rl%it,t^) = j{dq/2Tx)Rf\t,t^). The 
small and large r scaling forms of Eq. ( 148|) are 



r < : Rf{t,t^) ^ h^q^—^e--2^'\ (49) 

r > : 4'^)(t,tw) ^ l-t^q^^J^e--^^"'. (50) 

In analogy with the results for correlations, the exponents in the Gaussians differs by a 
factor of I between Eq. ( H9l) and Eq. (!50|) . For the construction of FD plots we need the 
step response function associated with Rq^\t,t^). This is obtained straightforwardly 
from Eq. (gH]), 

Xt\t,t^) = f dt'Rf\t,t') ~ 4t [l . (51) 
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It remains to work out the scaling of the symmetric part R'^\t, t„) of the response 
function. To calculate the Fourier transform of Eq. ( 120|) we use the identities Eq. (lA.Sp 
and Eq. (lA.9p . which produce 



i?f)(t,tJ~-e-2*cos(|'' 



(52) 



again with A = ^Jt'^ cos(g/2)2 + sin(g/2)2. We now perform the usual aging 
expansion — * oo with rjt^ and t^q^ fixed. But for Eq. fl52|) some care is needed: 
the leading terms in the expansion cancel. Therefore the first subdominant terms in the 
modified Bessel functions have to be retained. In this way one obtains the scaling 

^ (l - i^^^^) e-i»-^'"'«'. (53) 



The small and large r scaling forms of Eq. (1531) are 

r « : Rfit^t^) ^ -^^^^ (l - (54) 

r » t_ Rf{t.t^) ^ -^^T^^e'^^"'- (55) 

Once again the coefficient in the exponents of the Gaussians crosses over from ^ to ^ as 
we go from the regime r ^t^ to r t^. 

Integration of the expression Eq. fl53l) yields a simple result for the step response 
function associated with R^q\t,ty^), 

X^HMw) = f^dt'Rfit^t') ~ -^^^e-i(^-wA^)*.^ (56) 

Putting together the results from above yields the FD plots for the observables Uq 
in the FA model. In contrast to the case of the local observable Ui non-trivial limit plots 
do exist as long as we keep tq^ or equivalently \/tq constant as we increase t. We use 
normalized plots showing Xq versus ACg = 1 — Cq where 

If parameterized with t„, the slope in the plots directly gives the FDR Xq. For our 
wavevector observables we use Cq(t,t), Eq. (l43ll . to normalize the plots. The scaling 
expression for Cq{t,t^) is stated in Eq. (jlOl) and Xqit^tw) = xt\t^t^) + Xq\t,ty,) is 
given by Eq. flSTj) and Eq. fl56|) . respectively. Numerical evaluation of these quantities 
produces the FD plots shown in Fig. O 

We first note that observables Uq with \^q ^ 1 produce FD plots close to the 
equilibrium line Xq = ^ ~ Cq, similarly as is the case for the local defect observable nj. 
Very small values of q, \ftq <^ 1, on the other hand, yield essentially the same FD plot 
as the energy (g = curve in Fig. [2]). For intermediate wavevectors q the FD plots 
interpolate between these two extremes. The wavevectors q have an associated intrinsic 
lengthscale I = 1/q. Thus the observables Uq allow one to probe violations of FDT on 
any given lengthscale. While the coarsening dynamics equilibrate short length scales 
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Figure 2. Normalized FD plots for Fourier mode observables Ug in the d — 1 FA 
model. Symbols are simulation data at T = 0.08 in a system of linear size L = 700, at 
final time t = 4.87 x fO"^ and wavevector q = (27r/L)j with j = 0, 3, 6, 9, 12, 15, 18, 
21 and 27 (from bottom to top); note that t is given in microscopic units rather than 
units of 1/c, as used in the analysis. The dashed line represents equilibrium FDT. Full 
lines show the theoretical scaling predictions. The lowest curve is for q ~ and gives 
the FD plot for the energy. It is slightly curved, with slope increasing (in modulus) 
from -(1 + %/2) at the origin to X^, Eq. ([58]), on the right. 



i ^ ^/t, large length scales i ^ \/i remain far from equilibrium with a well defined 
violation of FDT corresponding to the case q = 0. 

The FD plots in Fig. [2] are unusual in two ways: first, they are non-monotonic and 
second, the step-response functions Xq can actually become negative. We remark that 
the non-monotonicity of the plots is genuine and not due to incorrect parameterization 
with t instead of t^ |112j . As for the local and global observables discussed above, 
the response function Rq comprises an asymmetric and symmetric contribution. In 

fg\ /r-7—i 

the regime r <^ t„, > dominates over Rq' < 0, see Eqs. ( H91 IMl) . One then 
verifies from Eqs. dH |49]) that Xq ~ Ri''\t,t„)/dt^Cq{t,t^) ~ 1. Equilibrium FDT is 
satisfied for t -^t^ and all FD plots have slope Xq ~ 1 close to the origin. This applies 
with the exception of g = itself; numerically, already for small but nonzero q the 
initial equilibrium regime is difficult to discern from the plots. In the opposite regime 
T ^ on the other hand, R^^ < dominates over Ri""^ > 0, c.f. Eqs. ([501155]), and the 
response Rq becomes negative. This causes the non-monotonic FD plots and negative 
step-response functions Xq- The resulting asymptotic FDR X^ for r ^ is given by 
Xq ~ Rq^\t,t„) /dt^Cq{t,t„) , and from the scalings Eqs. fH6 l [55]) . takes the value 

= -TT^^ ^ -3.307. (58) 

^ OTT — lb 

All FD plots in Fig. [2] thus end with the same slope X^ on the right, regardless of q. 
This is clearly observable in the plots with \/tq -C 1 but becomes difficult to detect when 
\/iq 3> 1. This reinforces the general point that observables which probe directly large. 
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non-equilibrated length scales are better suited for measurement of |125[ 1136] . 

Negativity of the response function Rg can be traced back to the activated nature 
of the dynamics. Perturbations coupling to the Ui affect the energy barriers for creation 
of defects on these sites. This effect and the associated response functions can be 

(s) (a.) 

decomposed into symmetric Rq and asymmetric Rg contributions, with the former 
trapping diffusing excitations but the latter accelerating their dynamics. In the regime 
T ^ the acceleration effect dominates, and the faster relaxation towards n{t) —>■ 
decreases the expectation values {ni(t)), thus the negative sign of the response functions 
Rq. We will see below that similar arguments apply to FA models in dimensions d > 1 
and also to the East model. 

The asymptotic scaling results derived in this section are confronted with numerical 
data in Fig. [2l The simulations are based on the continuous time BKL algorithm 
[138] . Step response functions Xg with g 7^ are measured using the no-field method 
of Chatelain [122[ 11231 1124] . while a specialised no- field method is employed for the 
q = case, see Sec. 12.41 below. These no-field methods enable us to construct properly 
parameterised FD plots with t„ G [0,t] as the running parameter. It is obvious from 
Fig.[2]that excellent agreement is obtained between the scaling results and the simulation 
data. 

2.3. Field theoretical approach 

In this subsection we describe how to calculate aging properties of the FA model by 
means of an effective field theory built from the Doi-Peliti |139i 1140] formalism for 
systems with stochastic dynamics. This approach to treat the FA model was first 
introduced in [21] and developed further in [23l[25. A summary of the results below, in 
so far as they relate to aging dynamics for dimensions d > d^ = 2, can be found in [125] . 
The aim of this section is to provide a detailed discussion of the calculations leading to 
those results, and to consider their extension to d < dc where a comparison with the 
exact results from the previous section can be performed. 

2.3.1. Effective theory and dynamical action We consider the "bosonic" version of the 
FA model [211 ISSl [23]: we allow for each site of the lattice to contain any number 
of excitations; that is, Ui = 0, 1, . . . indicates the occupation of site i. Since we are 
interested in the regime of low temperatures where the equilibrium density of excitations 
is low, removing the hard-core restriction on site occupation does not change the 
behaviour of the model significantly [211 ESI [21] . 

The derivation of the field theory for this version of the FA model was described 
in detail in Ref. [23] (see also [23] )■ The effective theory is defined for excitations on a 
lattice of sites in d dimensions with integer occupancy per site. As explained above, 
excitations can effectively diffuse between neighbouring sites with a rate proportional 
to c; neighbouring excitations can coalesce with a rate likewise proportional to c 
[2TI 11251 [23] • Excitations can also branch into pairs, but the rate for this process 
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goes as c^, and we disregard it in our analysis of aging at low temperatures. Using the 
Doi-Peliti formalism |139i I140[ I14H 1142] we can represent this dynamics in terms of 
complex fields 0(r,t), 0(r,t) with an action pTl [231 IM 1125] 

S = [ 0(9t-DV2)0 + A0(l + 0)02-no05(t), (59) 

Jr,t 

where D,X cx c. This is the field theory for the A + A A react ion- diffusion process 
[1411 1142] . In this representation the local occupation number operator n{r,t) is given 
by 

n(r,t)= [l + 0(r,t)]0(r,t). (60) 

The last term in ( 159|) indicates a Poisson distribution of random initial conditions of 
density uq. Correlation functions are calculated through the path integral: {A{t)) = 
J D(l)D(f)A[n{t)]e-^. 

2.3.2. Tree-level density and two-time correlations For d > = 2 the field theory is 
finite and a tree level calculation will give the correct behaviour at long times. Such 
a tree-level analysis amounts to solving the Euler-Lagrange equations. From (|59|1 we 
have: 

A *:^' 

{dt - DW^)<p + A(l + 20)02 - no5(t) = 0, (61) 



-(9t + DV')0 + 2A0(l + 0)0 = 0. (62) 



The average density is given by: 



n{t) ^ / (n(r,t)) = / (0(r,t)) = K-^(0o(t)), (63) 

Jr,t Jr,t 

where 0q(t) is the Fourier transform of (j){r,t), (f)q{t) = e~*'^'^0(r, t) [similarly for 
0q(t)], and we have used that averages with a factor on the left (i.e. at the latest 
time) vanish. Taking the expectation value of Eq. (1611) allows us to calculate n{t) at 
tree level: 

. ^ + A^-(^.>^ . ^ nm ^ ^ « 1. (64, 

From (l62l) we can get the tree-level propagator Gq{t,t^) = y~^(0q(t)0_q(tw))- To 
0(02), Eq. (|62|) can be integrated to give: 

0.(t„) = ^e-^«^(*-*»)0.(t). (65) 

Making use of the identity lim^-^o V~^{(j)q(t + e)(f)^q{t)) = 1, we obtain: 

G,(t,U = ^e-^«^(*-*-). (66) 
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Equations flMl) and (1661) now allow to compute two-time density correlations, 

C,{t,t„)= I e"'-^[(n(r,t)n(0,tw))-^(t)n(t„)] (67) 

J r 

= y-i(0,(t)0_q(t„) + 0,(t)0-q(tw)l^"Vo(tw)) (68) 
= ^"'(0q(t)0-q(tw)) + G,{t, t„)n(t„). (69) 
We now use f l6T|) to obtain the correlations at tree-level from {(f)qit)SS/S(j)q{t^)) = 0: 

dt„Cg{t, t^) + Xn{t^) {z + 2)C,(t, t„) - Xn\t^)G,{t, t^){2z + 1) = 0, (70) 

where 



Dq 



2 



DqH^. (71) 



An(t„) 

The general form of Cq{t,t„) that results is 

C,it,t„) = f{z)n{t„)G,it,t„), (72) 
where /(2;) obeys zf\z) + (2z + ^)f{z) — {2z + 1) = 0. The solution to this equation is: 

with z given by Eq. fITT]) . 

2.3.3. Tree-level response and FDR Consider now a perturbation /ig of wavevector g, 
cf. Eq. flTTI) . The corresponding change in the rate constants is Dq D6go + D(3hq 
and \q — > \5qQ + X[3hq. The change in the average density, to linear order in hq, is 
given by the linear order correction {4>q{t))^^^ to the expectation value of 0. This can 
be calculated at tree level by expanding {SS/Scj)) = to linear order in 5hq. 

' Dq^ 
_\n{t) 

Here we have used the fact that when the diffusion constant is not uniform, as is the 
case with a perturbation like the one here, the diffusion term in the action reads 
~ L,t^^ ■ {DV(t)-(tND) [125]. The response, Rq{t,t^), is given by Rq{t,t„) = 
T6{(f)q{t))^^y 5hq{t^), and obeys the differential equation: 



= dt{cf>q{t))^'^ + [Dq^ + 2Xn{t)] {4>q{t))^'^ - (3hq{t)Xn\t) 



- 1 



(74) 



^ 0, (75) 



dtRq{t, tw) + [Dq^ + 2Xn{t)] Rq{t, t^) - 6{t - t^)Xn\t) , ^^^^^ 

which implies the initial condition Rq(t^,tvi) = {z — 1) Xn"^ (t^) . By integrating the 
above equation we get: 

Rq{t, t„) = (^ - 1) An2(t)e-^^'(*-*-). (76) 

From equations (1721) and (!76l) we obtain the FDR for large t and t„ and for all 
momenta, 

X(ff)- _i-3 + 12z (.«!) 

"^ ' "^ {l + d.)zf{z) ] 1-1/z (^»1) ' 



Non- equilibrium dynamics of spin facilitated models 



25 



-2 




0.4 0.6 



Figure 3. Normalized FD plots for Fourier mode observables rig in the d — 3 FA 
model. Symbols are simulation data at T = 0.1 in a system of linear size L — 32, for 
final time t = 2 x 10^ and wavevectors q — tt/x with x = 1, 2, 2.4, 3, 3.2, 4, 5.33, 
6, 8, 12, 16 and oo (top to bottom). Full lines show the field theoretical predictions, 
accumulating on the equilibrium FDT line at large q. At the other extreme, the lowest 
curve coincides with the FD plot for the energy at q — 0, which is a straight line with 
slope —3. From Ref. [125j . Copyright American Physical Society. 



where z is given by Eq. ( 17Ti) as before. At any waiting time tw, for wavevectors q 
larger than 1/y/Dt^, FDT is recovered: Xg ^ 1. On the other hand, at small enough 
wavevectors, q -C 1/y/Dt^, the FDR becomes negative. In the limit g — > 0, we get the 
FDR for energy fluctuations, 

XE{t, tw) = Xg=o{t, t„) = -3. (78) 

From ( I77|) this is also the asymptotic FDR X^ = limt^^oo^i^t/t^^oo Xq(t,tw) for any 
wavevector. As for the one-dimensional case, a detailed comparison (see Fig. [3]) between 
these analytical results and direct numerical simulations of the FA model in dimension 
d = 3 shows very good agreement over the complete range of wavectors, from the 
straight line of equilibrium FDT for large q to the energy FD plot (straight line of slope 
Xe = —3) for small q. 

2.3.4- analysis and FDR for d < d^ For dimensions below the critical dimension 
dc the field theory needs to be renormalised. Our RG analysis here follows that of 
Ref. [1411 1142] . The theory only requires coupling constant renormalization. The 
Callan-Symanzik equations for two-time correlation and response functions, in terms 
of renormalized quantities, read: 

= [2{Dt)dDt + 2{Dt„)dDt„ - qdg + (3{gn)dg^ - dn^dna + A Cg{t, t„, uq, g^i, k), 
= [2{Dt)dDt + 2{Dt„)dDt„ - qdg + /3{gR)dg^ - duodno + id + 2)] Rg{t, t^; no, 5'R, k.), 
where we have made explicit the dependence of the functions on the initial density, 
the renormalized coupling constant, gn, and the arbitrary scale k which relates the 
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dimensional coupling A to the adimensional one g = {X / D) k,'^~'^ . The function 

is the exact beta- function [1411 I142j . These equations are solved by the method of 

characteristics (Till 1142] 



where (7R(t) is the running coupling |14H I142 J. In the asymptotic limit t ^ oo this 
approaches its fixed-point value, gR(t) — >• (= 27r for d = 1). 

We can now use the tree-level results ( !72l [76!) and Eqs. (!79| [80]) to obtain the 
renormalized tree-level correlation and response. To leading order in g^^ and for long 
times we make the replacement A — > Dg^K,'^~'^. Using (1791 ISUI) and setting c? = 1 we 
obtain for large t and t^: 

C,{t, t„) ^ / {DqH^) -1= f ^"j e-^«^(*-*»), (81) 



9k 



Dt \ t 



R,it, t„) ^ (DgX - 1) ^77=^ e-^«^(*-*-). (82) 

The time dependence of these correlations and responses is different from that for d > dc, 
see Eqs. fl72ti76p . The FDR, however, is the same as that given in fl77j) . The above 
expressions fail to reproduce the exact results obtained above for d = 1 (see Sec. I2.2.3p . 
Eqs. ( IHTl 1821) are the renormalized tree-level functions. In principle one could calculate 
loop corrections which will change the dependence on times and wavevector. It is likely, 
however, that the discrepancy with the exact results will persist. The problem is similar 
to what occurs with one-time functions [1411 1142] : the RG analysis is based on a d — dc 
expansion which is not quantitatively accurate at = 1. 

Simulations performed in dimensions d = 1 to d = 4 fully confirm our results, 
see Fig. |H While the d = 1 FD plot for the energy is curved and follows perfectly 
the theoretical predictions obtained in Sec. 12.2.21 the plots in higher dimensions are all 
compatible with the straight line of slope —3 calculated in this section. 



2.4- Simulation methodology 

Chatelain proposed a very useful no-field method for measuring response functions in 
spin systems with Markovian dynamics [1221 1123[ 1124] . The response is re-expressed as 
a multi-point correlation function, which has the advantage that it can be calculated 
from unperturbed trajectories. However, the efficiency of this method is poor for long- 
ranged perturbations: the statistics can become very noisy, especially for the collective 
quantities we want to study here. It turns out, however, that the response of KCMs to 
energy/temperature perturbations has special properties that make it possible to derive 
a more efficient no-field method. Specifically, in any KCM and any dimension [125] , the 
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l-CE{t,U)/CE{t,t) 

Figure 4. Normalized energy FD plots for the FA model in dimensions d = 1 to 4. 
Symbols show simulation data, all taken in the aging regime for T ~ 0.2 to 0.06 and 
t = 10^ to 10^. Full curves show a straight line with slope —3 and the exact limit plot 
in d = 1 with asymptotic slopes shown as dashed lines. From Ref. |125| . Copyright 
American Physical Society. 



step response can be written as 

XE{t, tw) = ^ [(1 - 2c) r dMt) + CEit, t) - CeH, t„) + Cyit, t) - Cyit, t^)] . (83) 

Here Cy(t,tw) = {l/N){{E(t)Y(t^)) - {E{t)){Y{t^))) is the normalized connected 
correlation between the energy E and the random variable 

Y{t)= / dt'U{t') with U = y^fi{n){ni-c). (84) 

i=i 

In Monte Carlo simulations the integrand in Eq. (18^ is constant between successive 
spin-flips and hence the integral becomes a sum Y{t) = J2ti<ti^i ~ ti~i)U{ti), with ti 
the sequence of updating events. Equation ( 18^ is easily implemented, especially in 
event driven code that keeps track of mobile sites i with fi{n) ^ 0. By recording 
histories of the energy Eii) and the observable l^(t), step response functions XE(t,t^) 
follow immediately from Eq. (!83|) . Importantly, one can obtain values of XE(t,t„) for 
arbitrary combinations of t, t^ by averaging over a stored collection of histories. This 
is crucial for the efficient generation of FD plots with fixed t and running t„. Without 
this method, it would have been impossible to obtain the collection of data shown in 
Fig.H 

In the remainder of this section we present a derivation of the no-field method 
Eq. (J83l) . It is convenient to use the operator formalism introduced earlier. We first 
consider the disconnected correlation DE{t,t^) = {E(t)E(t^)) which can be expressed 
in the form DE{t,t^) = {e\Ee'^'^ Ee^*""\po) , with |e) = J2n\'^) before and \pq) the 
initial state. Differentiating this with respect to t„ at fixed t (note that r = t — tw) gives 

dtMt,U) = {e\Ee'^^[E,W]e'^'-\po), (85) 
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where [A, B] = AB — BA denotes the commutator. For our derivation we will also 
need the response GE{t,t^) = T[S{E(t)) /6h(t^)] of the energy to the perturbation 
SE = —hY^-Ui. In the operator formalism this is 

G^(t,M = (e|Ee^-re^*"M. (86) 

The perturbation operator is = TdhW{h)\h=o where W{h) denotes the master 
operator in the presence of the field h. The basic idea for deriving a no-field method 
for the measurement of GE{t,t^) is simple: we try to form a linear combination 
aGE{t,t„) + bdt^DE{t,t„) = {e\Ee^^Xe^'-\po) with X = aV + b[E,W] so that off- 
diagonal components in X match the master operator W (there exist no coefficients 
a,b such that X is purely diagonal). Then this quantity can be expressed as a time 
derivative. It turns out that the relevant linear combination is a = 2 and b = —1. To 
see this let us work out the explicit form of V and [E, W]. From the master operator 
Eq. ([9]) we see that calculation of V only requires Tdhc{h)\h=o = c(l — c); the Glauber 
rate for activation in the presence of the field is c{h) = 1/[1 + e^^^~'^^]. So, 

N 

V = c{l- c) - - 2^*)- (87) 

The calculation of the commutator [E', VT], on the other hand, essentially reduces to 
[nj, Fj\ and this is [rij, Fj] = 5jjFj(l — 2nj), whence 

[E,iy] = ^F,/,(c-n,). (88) 

Based on the last two equations one verifies that X = 2V— [E, W] = {l—2c)W+U, with 
U = fi{ni — c) the (diagonal) operator corresponding to the observable U, Eq. (18^ . 
This is now a no-field method because 

2GE{t,t^)-dt^DE{t,t^) = (e|Ee^^Xe^*-bo) 

= (1 - 2c){e\Ee^^We^'-\po) + {e\Ee^^Ue^'-\po) 

= (1 - 2c)dt{E{t)) + {E{t)U{t^)). (89) 

Note that this equation applies to any KCM since we have not made any assumptions 
on the particular form of the kinetic constraint fi^n). We have only used that the 
unconstrained flip rates Wi{n) are Glauber rates for E = ^^Ui (a similar result can 
be derived for Metropolis rates). It remains to switch from disconnected to connected 
correlations in Eq. (l89|l and to multiply by so that all quantities are intensive. This 
gives 

2RE{t,t^) = (1 - 2c)dtn{t) + dt„CE{t,t^) + l((E(t)[/(t„)) + {E{t))dtAE{t^))). (90) 

Obviously RE{t,t^) = {l/N)GE{t,t^), CEit^t^) = {1/ N){DE{t,t^) - {E{t)){E{t^))) 
and n{t) = {l/N){E{t)). But the last term in Eq. (!90|) is in fact the connected correlation 
between E{t) and U{t^) since 

dtAE{t^)) = (e|We^*-bo) = -{e\Ue'^'-\po) = -{U{t^)). (91) 
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The identity {e\EW 



{e\U is most easily seen from Eq. fl88l) by multiplying both 



sides with (e| and using conservation of probability {e\W = a.s well as completeness of 
the projection state {e\Fi = (e|. Equation ( l9Ti) makes intuitive sense: according to its 
definition, Eq. (18^ . U is a. measure for how much the concentration over mobile spins 
(with fi{n) = 1) differs from the equilibrium concentration c. Thus, if {U{t^)) > 
there is an excess of mobile rij = 1 spins and {E{t^)) will decrease, and vice versa. 
Our no-field method Eq. fl83|) . finally, follows by integrating Eq. (|90|1 to obtain the step 
response XE{t,t^) = J^^dt'RE{t,t'). 

For completeness we add that this no-field method of course reduces to the standard 
FDT in equilibrium. This is most easily seen from Eq. ( |89l) . In equilibrium the 
term dt{E(t)) = drops out since the energy is stationary. Further, dt^DEit,t^) = 



-{e\EWe^^'-'-^E\p,^) = {e\Ue^^'-'-^E\p,^) = {U{t)E{t^)) = (E(t)f/(t„)), where 



the last equality expresses time-reversal symmetry. Thus Eq. ( l89l) retrieves the FDT 
GE{t-,t^) = dt^DE{t,t^) in equilibrium. 

3. East model 

In the second part of this paper we move to the East model, a directed variant of the 
FA model. It is defined on a one-dimensional lattice, i.e. a chain of spins, and differs 
from the FA model only in the choice of facilitation factor: 



Compared to = ni^i + rij+i in the FA case, this allows facilitation only from the 
left; the spin to the "East" of each defect is mobile. This seemingly harmless change 
has profound consequences. In particular, there are no finite groups of defects that can 
diffuse unimpeded over a defect-free (all down-spin) background. This is clear because 
the leftmost up-spin of any such group will never be able to change its state, due to the 
absence of a facilitating neighbour on the left. The East model therefore has much more 
cooperative dynamics, and its relaxation times correspondingly diverge more quickly as 
c decreases towards zero. In fact, when expressed in terms of temperature, relaxation 
times to equilibrium after a quench diverge as an exponential of inverse temperature 
squared, Intrei oc oc ln^(l/c) [271 [28 l 0.161 11171 1143] . compared to the Arrhenius 

dependence Intj-ei oc ln(l/c) ~ 1/T for the FA case. 

We consider again the dynamics after a quench from equilibrium at infinite 
temperature, where every spin is up independently with probability 1/2, to some low 
final temperature T. On timescales of 0(1) the system will reach a state where all 
up-spins are isolated. Further reduction in the defect density proceeds via thermally 
activated fluctuations, but now in a much more non-local manner than for the FA model: 
each up-spin can flip down only when, starting from the nearest up-spin on the left, a 
fluctuating "front" of up-spins has extended sufficiently far to the right to mobilize it. 
To be specific, consider an up-spin at site rf, and let the nearest up-spin on the left be at 
site 0. To relax this "domain" of length d, a front of defects needs to propagate from uq 



fiin) = 



(92) 




Figure 5. Domain size distribution in tlic East model in the first few plateaus of the 
dynamics. Open symbols are the theoretical prediction for the irreversible coarsening 
regime, filled symbols are from simulations. Circles: k = —1 (initial condition), 
squares: k = 0, diamonds: fc = 1, upward triangles: k = 2, downward triangles: 
fc = 3. The inset shows d^inPk{d) versus x = d/dmin to show the approach to the 
predicted scaling form P{x) (solid line) for large k. 



SO that it creates a defect at site d — 1, thus mobihzing n^. Naively one might suspect 
that the creation of this front takes an energy d — 1, with spins ni, n2, ■ ■ ■ , Ud-i simply 
flipped up in sequence. However, further reflection P, 1117] shows that the most 
efficient way of creating a front is hierarchical on lengthscales increasing as powers of 
two: after flipping up rii and n2, rii can be flipped down; then and n4 are flipped 
up, ris is flipped down, and finally 77,2 is eliminated by reversing the initial part of the 
process. Pictorially, one has 

10000 ^ 11000 ^ 11100 10100 10110 ^ 10111 10101 

11101 11001 10001 ... 

Effectively one creates a "stepping stone" at site 2, steps from there to site 4, and 
removes the stepping stone at site 2; this process can then be continued by creating a 
stepping stone at site 8 from the one at site 4 and so on. Keeping track of the maximum 
number of extra defects that exist at any one time, one shows that the energy barrier 
for flipping down the up-spin at site d (via creation of an up-spin at site d — 1) is 

k for 2^-^ <d<2\ (93) 

where k = 0,1,2,... |116[ 1117] . Such relaxation processes therefore take place on 
timescales t ~ exp(A;/T) ^ c~^. Importantly, in the limit of small c the ratio of any two 
of these timescales diverges, and the dynamics separates naturally into "stages" labelled 
by the value of k. We can include in this picture the initial relaxation after the quench, 
where domains of length d = 1 disappear without activation; this is stage k = 0. 

To understand the non-equilibrium dynamics in more detail, it is useful to view it 
as a coarsening process: when the up-spin at the right boundary of a domain disappears. 
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Figure 6. Defect density n{t) versus vt for tliree different temperatures. Symbols 
represent the simulation results. The theoretical prediction pi8p is shown at finite T 
with dotted lines, while its T — + limit is shown as a full line. 



this domain merges with the neighbouring domain on the right. Because the rate for 
disappearance of an up-spin depends only on the length d of the domain it bounds, 
it is clear that no correlations between neighbouring domain lengths can ever develop 
during such a process |116[I117] . For low c, the evolution of the domain size distribution, 
P{d, t), can therefore be calculated without approximation by an "independent interval" 
approach. In stage k of the dynamics one finds the equation of motion (for d > 2^) 



^P{d,t)= J2 Pid-d',t) 



dt 

2'=-^<d'<2'' 



(94) 



The sum runs over the "active" domain sizes d' that are eliminated during stage 
k. The rate for creating longer "inactive" domains is then the product of the rate 
[—{d/dt)P{d',t)] at which domains of size d' disappear, and the probability of having 



a domain of length d — d' on the right. As shown in Appendix B , this equation can 
be solved to relate the domain size distributions at the beginning and the end of 
stage k, and this relation can then be iterated numerically, starting from the initial 
P{d, t = 0) = 2~'^ |116j . The results are shown in Fig.[5]for the domain size distributions 
Pk{d) in the first few "plateaus" of the dynamics. Here we label each plateau by the 
stage which it follows, i.e. plateau refers to times 1 <^ t <^ c~^. More generally we 
will use the notation 

z/ = Tlnt, A;=[z/J, a = u — k, (95) 

so k denotes the plateau that has been reached (and the stage of the dynamics that has 
finished) at time t. To avoid confusion we note that plateaus were labelled by + 1 in 
previous work [ 1161 1117] ; in our convention, the initial condition corresponds to plateau 
k = -1. 
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As Fig. O shows, the predicted domain size distributions agree well with simulation 
data at low T. From the domain size distribution one can predict in particular the 
defect density, giving n(t) = 1/dk where dk = '^d^^kid) is the average domain 
length in plateau k. Again this agrees well with simulation data (Fig. [6]), though the 
sharp steps that are predicted for T — >■ in a plot of n vs at the integer values 
u = 0,1,2, .. . are obviously rounded at nonzero T. By analysing the dynamics in more 



detail (see [Appendix B| and Eq. (11181) below) also this effect can be predicted, and then 



gives a very good description of the data. 

A further observation from Fig. [5] is that the domain size distributions approach 
a scaling form for large k. This is obtained by scaling with the minimum domain size 
present at time t, 

draUt) =2^ + 1, (96) 

so {t)Pk{d) is plotted versus x = d/djain{t)- The limiting distribution P{x) 

can in fact be worked out explicitly |116l 1117] : its Laplace transform is 1 — exp[— Ei(s)] 
where Ei(s) = dzexY){—zs) / z is the exponential integral. Expanding in powers of 
Ei(s) and inverting the transform gives then 

1=1 

where the functions fi{x) can be defined recursively by successive convolutions, 

f,{x) = 5{x), f,{x) = Q{x-l)-, fi^,{x) = {fi*fr){x). (98) 

X 

Here Q{x) and 5{x) are the standard Heaviside step and Dirac delta functions. Explicitly 
one has f2{x) = 0(x — 2)(2/x) ln(x — 1), while for / > 3 the fi{x) cannot be expressed 
in any simple closed form. Fortunately each fk{x) contains a factor Q{x — k) so for 
realistic ranges of x only the first few functions are needed. The mean value of x across 
the distribution (!97|) is exp(7) = 1.78..., with 7 being Euler's constant. Converting 
back to the unsealed defect density gives then 

n{t) = e-^2-^ (99) 

for large k. Because of the logarithmic dependence of on t, this scales as n{t) ~ t^-^'"^^: 
the East model exhibits anomalously slow coarsening, with the decay exponent Tin 2 
that governs the decay of the defect density decreasing towards zero at low temperature. 

We will refer to the above analysis of the East model dynamics after a quench to 
low c (or equivalently low T) as "irreversible coarsening" . It becomes exact in the limit 
c — > taken at fixed stage k, i.e. with typical domain sizes d kept fixed and of order 
unity. This way of taking the limit ensures that up-spins do indeed flip down irreversibly: 
the probability of observing, at some given point in time and within a domain of length 
d, an up-spin as part of a fluctuating up-spin front is 0{dc) and can be neglected in 
the limit. Also, even though within a given stage k there can and will in general be 
multiple relaxation timescales because active domains cover a range of lengths and can 
be eliminated by passing through different sequences of intermediate configurations, the 
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timescales for different stages do always separate by much more than this spread for low 
c. 

Both of the above simplifying properties no longer apply if we wish to look at 
typical domain sizes of the order of the equilibrium domain length d^q = 1/c. We 
would then take c small at fixed 6 = dc, rather than at fixed d. For 6 of order unity 
and larger, this gives information on the eventual crossover to equilibrium behaviour. 
For the purposes of this paper we will only be concerned with the earliest stages of 
this crossover, 5^1. This should connect appropriately with the large-c? limit of the 
irreversible coarsening regime, where from (!93|) the timescale for relaxation of a domain 
of length d is t ~ c~'' ~ ^-in(i/in2 ^i/(Tin2)_ simple scaling hypothesis, supported 
by numerical evaluation of first passage times for the elimination of domains of large d, 
then gives the following picture [117j . The splitting into discrete stages of the dynamics 
is lost because the spectrum of relaxation times within the stages becomes so broad 
that a clear separation no longer exists. However, there is a now a continuous form 
of timescale separation: the time t for eliminating a domain of scaled size 6 is again 
t ~ (^/c)^/^^'"^^-*. Comparing two different domain sizes 5i and 82 then shows that one 
always has t2/ti ^ 1 in the limit T — 0, as long as 62/61 > 1. As a consequence, each 
up-spin has a sharply (on the spatial scale (i ~ 1/c) defined "equilibration zone" to its 
right, of size 6 ~ ct~'^^^'^. As t increases these zones grow; when an equilibration zone 
reaches the next up-spin to the right the latter is eliminated and the two domains either 
side merge. Put differently, the dynamics consists of always eliminating the shortest 
domain and "pasting" all of its length on to the domain on the right. This "paste- all" 
model has been studied independently in the literature |144] . Remarkably, the domain 
size distribution it produces in the scaling limit is exactly identical |117j to the one we 
found in Eq. fl97j) . This shows that the irreversible coarsening regime of the East model 
connects smoothly to the paste-all regime just discussed, as it should: in the large-/c limit 
of irreversible coarsening the domain size distribution approaches ( 1971) . and this form is 
preserved - apart from a trivial scaling with the minimum domain length d^i^ - in the 
paste-all dynamics. The paste-all regime has two advantages for theoretical analysis: 
only the smallest domain length is "active" and hence the temporal order in which 
domains disappear is deterministic; this is not the case in the irreversible coarsening 
regime, where two neighbouring domains of different lengths can be simultaneously 
active so that either can disappear first. Secondly, the role of the clock is played by 
c^min, so that effectively diverging timescales, 

t^d^^-'\ (100) 

can be considered. Even where predictions cannot be extracted analytically, they are 
then relatively easy to obtain by numerical simulation. 

After this overview of the theoretical tools used to analyse the low-T out-of- 
equilibrium dynamics of the East model we turn as in the FA case to specific correlation 
and response function pairs, first for local observables, then for the global one (which 
is just the energy), and finally for non-local observables. Because we need to deal with 
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Figure 7. Local correlation function difference from simulations, for i^t as indicated 
and T = 0.15 (symbols). Dotted lines: leading order prediction (|103p . Solid lines: 
theory including quasi-equilibrium corrections (jll3p . Note the logarithmic y-axis; on 
a linear axis the correction terms would be difficult to discern. 



two-time functions in all cases it will be useful to refine our notation from ( l95i) to 

Ut = T\iat, kt=[ut\, at = ut-kt, (101) 

with analogous quantities defined for t^ (giving u^, k^, a„) and the time difference 
T = t — t^ (giving kr, a-r). Simulations are carried out for system sizes N = 250 
and for temperatures T = 0.2, 0.17, 0.15, with corresponding equilibrium up-spin 
concentrations of c = 6.69 x 10~^, 2.78 x 10"^ and 1.27 x 10^^. For numerical 
convenience we mainly restrict ourselves to waiting times t^ > 1 when simulating two- 
time quantities. 

3.1. Local correlation and response 

We begin with the local correlation function Co(t,t^) = {ni(t„)ni{t)) — n{t^)n{t). 
Within the irreversible coarsening regime the time-dependence of this quantity is easy 
to determine: any spin that is up at time t must also have been up at time t^ since to 
leading order in c spins never flip back up. This implies {ni(t^)ni(t)) = {uiit)) = n(t) 
and so 

Coit,t^)=n{t)[l-nit„)]. (102) 

Simulation data are compared with this prediction - which is, trivially, exact for t^ = t 
- in Fig. [71 To make our later understanding of the FD behaviour easier we plot the data 
in the same form as required in an FD plot, showing the correlation function difference 
ACo(t,tw) = Co{t,t) — Co{t,tv/) for which the prediction (11021) takes the form 

ACo(t,tw) = n{t)[n{t^) - n{t)]. (103) 
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Figure 8. Local susceptibility xo{t, tw) in the quasi-equilibrium regime, scaled by cn(t) 
and plotted against Ur = Tln{t — t„). The steps show the theoretical prediction (|108p 
for c ^ 0. Squares correspond to T = 0.15, circles to T = 0.2. For both temperatures 
we show data for vt = 0.25, 0.5, 0.75, 1.25, 1.5, 1.75, 2.25, 2.5 and 2.75 (from left to 
right). 



This is plotted against t„ or, more precisely, z/^. The theoretical estimate fllU3l) 
rationalizes the simulation data quite well, but there clearly are small corrections for 
nonzero c. These are most noticeable when t and t^ are close together inside the 
same plateau of the dynamics; Eq. fll03p then predicts a very small correlation function 
difference because n(t^) ~ n(t). The actual values are somewhat larger because there 
is an additional quasi-equilibrium contribution in this regime. This correction can be 
most easily determined via FDT from the local susceptibility Xo? to which we turn next. 

To understand the effect of a field hi applied at time t^ on the evolution of spin 
Hi, let us take the history ni^i{t') of its facilitating left neighbour as given in the time 
interval t' = t^ . . .t. Spin Ui then behaves like an isolated spin except for the fact 
that it cannot change state whenever nj„i(t') = 0: the time interval available for it to 
relax from the state ni{t^) is reduced from t — t^ to dt' ni^i{t') . Since with Glauber 
rates (EI) the relaxation time of an isolated spin is unity, we thus have 

{uiit)) = c' + K(t„) - c'] exp (- [ dt'n,.i{t')] . (104) 



Here, c' is the perturbed up-spin density at site i. To linear order in the field strength 
hi it can be expanded as 

^ c + l3hic{l-c) + ... (105) 



1 + e'3(i-^') 

Inserting into (11041) and differentiating w.r.t. phi, the initial value Uiit^) drops out. 
After averaging over the history of nj_i one thus gets for the local susceptibility the 
general expression 



Xoit,t^) = c(l - c) 



(exp - / dt'm.iit') ) 



= cil-c)[l-r{t,t^)]. (106) 
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For low c the factor 1 — c can be dropped, and it remains to understand the behaviour 
of the "relaxation integral" r{t, t^). The average over the dynamics in its definition can, 
for a large system, be replaced by an average over sites j = i — 1: 



Sites where a spin has been up for a total time larger than unity make a negligible 
contribution, so r(t, t„) is to a good approximation the fraction of sites with spins that 
have remained persistently down between t^ and t; in other words, Xo(t,tw)/c is the 
fraction of sites that have been up at some point between and t. From this we can 
predict its dependence out — t^ in the quasi-equilibrium regime, where t^ and t remain 
in the same plateau of the dynamics so that no domains disappear and the up-spins 
defining them can be regarded as fixed. As long as r = t — t„ ^ 1 but r ^ c~^, 
only these fixed up-spins contribute to Xo', since these also determine the overall up-spin 
density, one predicts Xo/c = n(t^). Once r increases past c~\ the right neighbour of 
each fixed up-spin will have had time to flip up once or more, so that Xo/c = 2n{t^). 
This process then proceeds: for r ^ c~^, the relation ( l93|l between energy barriers and 
distances implies that the 2^ — 1 spins to the right of each fixed up-spin will have been 
up at some point. (The "—1" arises because we are considering flipping spins up\ the 
energy barrier for flipping a spin at distance d down, as written in ( |93l) . is the same as 
that for flipping a spin at distance d — 1 up.) Adding the contribution from the fixed 
up-spin back on gives overall 



(Either n[t^) or n{t) can be used in the prefactor since by assumption we remain in 
the same plateau.) This quasi-equilibrium prediction is checked against simulations 
in Fig. [HI with good agreement. Our use of the term "quasi-equilibrium" is justified 
by the fact that exactly the same step-like structure in xo would be observed in 
equilibrium, except for the replacement of the density of "fixed" up-spins, n{t), by 
the appropriate equilibrium value c. The simple scaling fllOSp of the step heights with 
powers of two does not seem to have been noticed in previous analyses of the equilibrium 
dynamics p fTlSlfTTT] . 

The quasi-equilibrium regime lasts while t and t^ remain in the same plateau of the 
dynamics, i.e. kt = k^. In the limit of low c we have kt = max(fcw, kr), so an equivalent 
condition for quasi-equilibrium is kr < k^. When t is held fixed as in Fig. [HI the limit 
kr = k„ corresponds to t ^ t^^ ^ t/2 and so, again for low c, kr = h. We thus expect 
to see kt steps in xo within the quasi-equihbrium regime, as confirmed by the data in 



Beyond the quasi-equilibrium regime we can continue to identify Xo/c with the 
fraction of the chain swept by up-spins, but this now becomes more difficult to calculate 
as domains and the corresponding "swept out" areas merge. Within the paste-all regime, 
however, it is a relatively simple matter to get accurate predictions by simulating the 
paste-all dynamics and keeping track of the swept out areas within each domain. The 




(107) 





Fig. [a 
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Figure 9. Local susceptibility xo in units of c, plotted against 9 — n{t^)/n[t) to show 
the out-of-equilibrium regime. From bottom to top, we have vt = 0.5, 1.5 and 2.75. 
Horizontal lines indicate the predicted height of the quasi-equilibrium regime for the 
three cases. Solid line: theoretical prediction for the paste-all limit {vt ^ 1). 



clock variable for this simulation is (imin(t)/(imin(^w)- Since in the paste-all regime the 
up-spin density scales with the inverse of (imin, we can equivalently write 

Xo(t, tw)/c = T{e), d = n{t^)/n{t). (109) 



We explain in Appendix F how the initial part of the scaling function can be calculated 
analytically, with the result 

J^(0) = e"^ (^1 + In - ^ In^ 0^ , for < 2. (110) 

This is in very good agreement with the numerical values extracted from the paste- 
all simulations, as demonstrated in more detail in Fig. [12] below. We show in Fig. [9] 
the overall prediction for Xo/c versus n(t^)/n(t) and compare with numerical data. 
Even though the simulations are performed within the first few plateaus, i.e. far from 
the paste-all regime, the data for increasing ut clearly do approach the theoretical 
prediction. The initial discontinuity of the theory represents the contribution from 
the initial quasi-equilibrium regime, which in the limit of low c shrinks to the point 
n{t^)/n{t) = 1. The value of Xo/c at this point is exp(— 7) ^ 0.56: each domain 
contains a swept-out equilibration zone of length d^i^{t^), and the density of domains 
is n{t„) = exp(— 7)/(imin(t„), giving a swept out fraction d„iin{t„)n{t^) = exp(— 7) of 
the chain. This is entirely consistent with extrapolating to the same point from the 
quasi-equilibrium regime: as discussed above, the latter ends where kr = k^- At this 
point, from (11081) . Xo/c = n(t^)2'''~ which for large k^ approaches exp(— 7) from ( l99l) . 

Before combining correlation and susceptibility data to get the FD plots we return 
to the quasi-equilibrium corrections to the local correlation function. To determine these 
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we use the following relation between local susceptibility and correlation: 

Xo{t, tw) = n{t)[l - n{t^)] - Co{t, t„) + c[n(t„) - n{t)]. (Ill) 



This is derived in Appendix C remarkably, it is exact and applies not just to the East 
model but in fact to all spin models with directed constraints. We used it to obtain 
the numerical data for the susceptibility in the simulations, but also checked it against 
direct susceptibility measurements. From the theoretical point of view we note that the 
first two terms on the r.h.s. of fillip cancel to leading order because of fll02p : this of 
course must be so because from (11061) the susceptibility never becomes larger than c. 
Because the 0(c)-corrections to the correlation function are difficult to estimate directly, 
Eq. (imp would not be useful to deduce accurate estimates for xo- We can, however, 
turn it around to get an expression for the correlation function (difference): 

ACo{t,t^)=n{t)[n{t^)-n{t)]+Xo{t,t^)-c[n{t^)-n{t)]. (112) 

The last two contributions on the r.h.s. are the exact corrections to the irreversible 
coarsening estimate (llOSp . They will be significant only in the quasi-equilibrium regime 
n{t^) ~ n(t), where furthermore the last term is negligible compared to the second one. 
The remaining dominant contribution to ACo(t, t^) is simply Xo(^5 ^w), as expected from 
FDT for dynamics in the quasi-equilibrium regime. Using (llOSp gives explicitly 

ACo(t,tw) ~ n{t)[n{t„) - n{t)] + cn{t)2^- . (113) 

This improved prediction is included in Fig. [7] above and now shows very good agreement 
with the numerical data. As expected, the correction term only affects the outcome when 
t and tw are close together and the correlation function is still close to its equal-time 
value; elsewhere it is negligible. 

We can now combine the above results to obtain the FD plot for the local observable 
in the East model. The raw results, with both axes scaled by Co(t,t) to produce 
the normalized two-time correlation Co(t,t^) = Co(t,t^)/Co{t,t) and susceptibility 
Xoit^tw) = x{'t^tv/)/Coit,t), are shown in Fig. [10, for fixed t with t^ varying along the 
curves. The data show clearly the initial quasi-equilibrium regime. The curves depart 
from this where tw becomes small enough to have left the plateau that time t is in. The 
kinks in this out-of-equilibrium part of the curves, where points accumulate, correspond 
to being well within a plateau; the next line segment begins when this plateau is being 
left and t^ moves (backwards) through the previous stage of the dynamics. The segments 
between the kinks appear straight to a good approximation, but close inspection of the 
data for larger Ut (as included in Fig. [9l and Fig. [H] below) does show a small amount of 
curvature, as observed also in the not dissimilar triangular plaquette model [126] . The 
slopes of the segments, which give the local FDR Xo{t,t„), are indicated in the plot. 
They are significantly smaller than unity but do not, in this raw form, follow a clear 
pattern. 

To determine a more appropriate scaling for the FD plot we recall from f llOOp that 
Xo/c is a scaling function of n{t„)/n{t), at least in the paste-all limit. For the correlation 
function difference we see from f ll03p that ACo(t, t^)/n'^(t) = n(t^)/n{t) — l is a function 
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Figure 10. Normalized FD plot for the local observable for T = 0.2 (full symbols) 
and T — 0.15 (open symbols). The numbers indicate the value of the fluctuation- 
dissipation ratio Xq in the various straight-line segments. For T one expects the 
plot to depend only on kt = [t'tj , while for nonzero T the sharp change at integer 
values of vt will be replaced by a smooth crossover. This is visible in the additional 
short linear segments just beyond the departure from quasi-equilibrium (dotted line) 
in the plots for T = 0.2 and ft = 0.75, 1.75. 



of the same scaling variable. (Since we are concerned with understanding the scaling 
in the non-equilibrium regime, we ignore here the small quasi-equilibrium corrections.) 
These considerations motivate one to plot, as we do in Fig. [11] (top), xo{t,t^)/c vs 
ACo{t,t^)/n'^{t). The curves for different temperatures now collapse reasonably well, 
except for small shifts caused by the presence of the initial quasi-equilibrium regime. 
If these are removed by subtracting from both xo and ACq the relevant contribution 
cn{t)2^\ as done in Fig. [11] (bottom), the collapse is really quite good. 

Looking more quantitatively at the local FDR and comparing with the FA case, 
the most striking feature is that for the East model Xq is never negative, even though 
the dynamics is activated. This is not just because any negative values are hard to see 
on the scale of the FD plot, as was the case for the FA model at late times. In fact, 
from fll06p one has 

Ro{t,t^) = —^Xo{t,t^) = c{l-c)(n,^i{t„)exp(^- rft' n,_i(t')^ ^ , (114) 

and because this local impulse response is always positive, the same is true of the local 
FDR. The magnitude of Xq, on the other hand, is very small: the scaling of Fig. [TT] 
shows that Xq is roughly of the order of c/n'^{t) for small c. To be precise, the expected 
scaling outside the trivial quasi-equilibrium regime is 

Xo{t,t^) = ^2(^) "^^t,!^^- (115) 

This contains the overall pref actor from the scaling of the FD plot. The remaining 
coefficient Sk^^u^ depends only on the plateau kt in which t is located, and the stage z/„ 
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Figure 11. (Top) Scaled local FD plot for T = 0.15 (open symbols) and T = 0.17 (full 
symbols) for different values of vt- (Bottom) Same with quasi-equilibrium correction 
subtracted; the data for i/t — 2.5, 2.75 are shown separately in the inset for better 
visibility. 



of the dynamics that tw is traversing. Note here that the FDR is defined only for integer 
z/„ in the hmit c ^ 0: for all other values, t^ lies within a plateau and Cq and A^o 
are both constant. Table [T] shows the values of S determined from our numerical data. 
They show a relatively weak dependence on temperature, consistent with the expected 
approach to nonzero limits for T — > 0; as temperature decreases, also the expected 
independence from the non-integer part at of ut becomes manifest. 

In our simulations we are necessarily restricted to exploring only the first few 
plateaus. Theoretically, we can also study the interesting paste-all limit where domain 
sizes are much larger than unity but still small compared to the average equilibrium 
domain length 1/c; this could be viewed as the true asymptotic non-equilibrium 
regime. We have included in Fig. [H] (top) the paste-all prediction for the scaled FD 
plot. It is obtained by combining the predicted scaling fllOQp with ACo(t,t^)/n'^{t) = 
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0.0257 
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Table 1. Numerical values of Skt,v^ for at — 0.5 (top) and at — 0.75 (bottom). 
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Figure 12. Paste-all prediction for xo{t,'tw)/c (dashed line) and [n^(tw)/c]Xo(i, iw) 
(solid line), plotted against 6 — n{tv,)/n{t). The lines show the results of numerical 
simulations of the paste-all dynamics; the circles give the closed-form prediction for 
9 <2. 



n{t^)/n{t) — 1 =6 — 1. Beyond the quasi-equilibrium regime the hmiting scaled FD plot 
is a smooth curve (although from the theoretical analysis in Appendix F one expects 
discontinuities in higher-order derivatives at integer values of 6). Within a mean-field 
picture, this would be interpreted as arising from an infinite hierarchy of relaxation 
timescales, each "responsible" for an infinitesimal segment of the FD plot, as found for 
instance in the infinite range Sherrington-Kirkpatrick spin glass model [62]. Remarkably, 
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Figure 13. Energy variance CE{t,t) versus vt for three different temperatures. The 
soHd line shows the theoretical prediction for the limit T ^ 0. 



even though the East model with its purely local and directed interactions is very far 
from mean field, this is exactly the scenario here. Along the x-axis of the plot we have 
9 — 1 = n(t^)/n(t) — 1 = dminit) / drainitv/) — l] two different points 9 and 9' > 9 then 
correspond to waiting times with, from fllOOl) . a ratio t'^/t^ = [9' /9Y^'^^^^'^^ that diverges 
exponentially for T ^ 0. 

We note finally that for the closely related triangular plaquette model, approximate 
relations for correlation and response functions have recently been derived by assuming 
that relaxations in the different stages of the dynamics are independent of each 
other |126] . The resulting factorization of the local correlation functions does hold in the 
East model for small c, as already discussed in [126] . The local FDR was also predicted 
to be independent of the later time t. To check the quality of this approximation at 
least in the paste-all regime, we differentiate ffTO w.r.t. /^CQ{t,t^)/n^{t) = 9 - lio 
get Xo(t,tw) = [c/n^{t)]J^'{n{t^) /n{t)). Independence of t would then require that 
[n^(t„)/c]Xo(t, tw) = [n{tw) I n{t)]^ J-'' {n{t^) / n{t)) be independent of n(t), i.e. constant. 
As Fig. [T2l shows, this is a reasonable approximation when n{t^)/n{t) is large, but for 
smaller values there are deviations which are in fact oscillatory in n{t^)/n{t). In the 
earlier, irreversible coarsening stages of the dynamics, t-independence is also not a very 
accurate approximation. To check this one can consider, instead of Sk^^u^ from Tabled! 
the combination Xq/c = Skt,u„n'^it). Fixing for example z/^ = and comparing kt = 
and kt = 1 one finds values which, while closer than those of 5*^^^,^^ itself, still differ by 
a factor around three. 

3.2. Global observable 

We now switch from the local observable Ui to its global analogue E = J^i^ij 

FA case, this is just the energy function of the system. The associated correlation and 

response are defined in ( |29l 1301) . 
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kt = 2 
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0.000580 


0.0188 


T = 0.17 


0.00376 
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0.000148 


0.0167 


r = 0.15 


0.00367 


0.00232 


-0.000110 


0.0160 


Theory 


0.00385 


0.00232 


-0.000233 


0.0153 



Table 2. Numerical values of energy correlation Csit, tw) by plateau of the dynamics, 
labelled by and kt > few Data for three temperatures are shown alongside the 
theoretical prediction for T — s- 0. In this limit, few = — 1 corresponds to iw = 0; 
the simulation data were in fact taken at a nonzero corresponding for each T to 
= —0.8, hence the slightly more noticeable deviation from theory in the first column. 



We start in Fig. [13] with simulation results for the equal-time correlation CE(t,t). 
Plotting against ut reveals that this has a plateau structure similar to the defect density 
n(t), and for c —>■ one expects sharp transitions to occur at integer values of z/^. This 
can be confirmed by a calculation within the irreversible coarsening regime, as outlined 
in Appendix D It produces values in very good agreement with the simulation data, 



as indicated graphically in Fig. [131 The theory shows that, not unexpectedly, energy 
fluctuations are related to the variance of the domain size distribution, via 

CEit,t) = ^^. (116) 

The same calculational approach can in fact also be used to predict the two- 
time energy correlations CE(t,t^), which have a plateau structure in both t and t^. 
Table [2] summarizes the predictions as a function of fc„ and kf. good convergence of 
the simulation results to the theory is observed as c is reduced. The only exception is 
(few = ^,kt = 2), where the predicted negative correlation is too small to be accurately 
determined from the simulation data. That such negative correlations could arise can 
be motivated by considering e.g. configurations with atypically many short domains at 
time t^, i.e. a high value of the energy; in the next stage of the dynamics these will have 
to merge and so can lead to unusually many large domains at a later time t. 

One may wonder whether negative energy-energy correlations are only a quirk of 
the first few stages of the dynamics. To check this we calculate in Appendix D the 
behaviour in the paste-all limit, where CEit,t^)/n(t) becomes a scaling function Q{9) 
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1 2 3 4 5 

Figure 14. CE{t,tyf)/n{t) from simulations, plotted against n{t^)/n{t). The various 
dashed lines show simulation data for vt — 0.5, 1.5, 2.5 and 3.5 (top to bottom). 
Results are displayed for temperatures T — 0.15 (except for vt — 3.5), 0.17 and 0.2 
but the variation with temperature is almost invisible on the scale of the graph. As vt 
increases, convergence towards the curve predicted in the paste-all regime (solid line) 
is observed. For finite vti i-e. in the irreversible coarsening regime, only the plateau 
values can be predicted theoretically (see Table [5]) ; these are shown by the symbols 
and agree closely with the simulation results. 



of = n{t^)/n{t). In terms of the functions fi{x) defined in fl98|) this scahng function 
reads 

G{e) = e-\e + + d9' fiie')^ . (ii7) 

It is plotted in Fig. [T3] and shows that (weak) negative correlations do persist even in 
the paste-all limit. The figure also shows graphically the simulation data, with symbols 
indicating the predicted plateau values from Table [21 In spite of the relatively low values 
of i>t reached, the data are clearly already moving towards the paste-all limit. 

The excursion of CE(t,t„) to negative values as illustrated in Fig. [H] implies in 
particular that the two-time energy correlations are not monotonic in t„ (at fixed t) as 
one would usually expect. An additional source of non-monotonicity is the behaviour 
around z/^ = 0, where in both simulations and theory we find that Ce starts to rise 
again as t^ decreases. We attribute this unusual behaviour to the relatively large value 
of the equal-time correlation CE(t^,t„) = 1/4 in the initial plateau (/c„ = —1). 

Next we look at the energy susceptibility XE(t, tw)- By definition this is the response 
of the defect density n{t) when the equilibrium up-spin concentration is changed at time 
from c to c', given in (11051) . To gain some qualitative insight we exploit the plateau 
structure of n{t) in the unperturbed dynamics. With successive relaxations taking place 
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Figure 15. — X£;(i,0), plotted against vt for tliree temperatures. Symbols: simulation 
data; lines: theoretical prediction, Eq. (I122p . 
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Figure 16. Negative energy susceptibility — xb(^:*w) for T = 0.15 and several vt, 
plotted against tv^/t. The proportionality to (1 — iw/i) expected from theory, Eq. (|122p . 
is well verified. 



on timescales 



In 



^, A; = 0, 1, 2, ... we can write (see Appendix B ) 



n{0) 



oo 

E 

k=0 



:ii8) 



where the functions gk{-) describe the relaxation within stage k and are independent 
of c to leading order. Each gk{C) increases from zero at C = and exponentially 
approaches a finite limit for large C,. For the first few stages one can calculate explicitly 
(see Appendix B ) 

1 



^?o(C) 



(1 - e-0 



(119) 
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Figure 17. The exponent a as obtained by fitting the c-dependence of XE{t,0) for 
fixed ft, compared with the theoretical prediction X-e(^)O) ~ c" with a — 1 — at. 



9,{0 = l{l-e-<'^), (120) 

S^(0-Yi^l-e-'"') + 'i{l-e-'")- (121) 

Now when the field is switched on, the effective time interval appearing in the argument 
of these functions will be replaced by cH^ + {c'Yit — t^). (For t„ 7^ this will not 
necessarily be exact except when the g^^ are single exponentials, but should give a 
reasonable approximation nonetheless.) Using ( 11051) and differentiating w.r.t. (3h then 
produces the following estimate for the energy susceptibility: 

XE{t, K) = -(1 - twA)n(t) kchg'kich). (122) 

k>0 

Each of the terms in the sum produces a "bump" in XE{t,0) at t ~ c"^ with a 
height of order unity; the dependence on t„ is only through the simple factor 1 — t^/t. 
Simulation results are in very good agreement with this prediction: see Fig. [T3 We plot 
~XE(t,0) because the susceptibility is predicted to be negative - recall that the gk{-) 
are increasing functions - in a clear signature of activated dynamics. Also the predicted 
linear dependence on is well verified by our data, as shown in Fig. [161 

The bumps correspond to the various stages of the dynamics; indeed, within the 
plateaus between these stages the defect density is independent of c for small c and so one 
expects the susceptibility to vanish for c — > 0. More precisely, if we write Tint = kt + at 
as usual with < < 1, we have the scaling 

-Xi?(t,0) -c-^, a = l-at. (123) 

To see this, note that all terms with k < kt in the sum (11221) are exponentially 
suppressed because the argument of gk{-) is large. In the remaining terms, on the 
other hand, the function argument vanishes for c — > 0; since the derivatives g'f.{0) 
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Figure 18. Top: Energy FD plots for i^t = 0.7 (left) and 1.7 (right) for tliree different 
temperatures. Bottom: Same plots with each axis rescaled by 1/c; fits to the linear 
parts give slopes of —2.55(2) and —3.77(5), respectively. 

are constants of order unity, each term scales as c^t. The dominant contribution 
thus comes from k = kt + 1, giving —XE{t,0) ~ ^h+i^-kt-at _ (d-^t^ claimed. 
Intuitively, the dominant contribution to —Xe is from the bump on whose uphill flank 
t is situated. The scaling f ll23p is checked against numerical fits of the c-dependence in 
Fig. [T71 Qualitatively we clearly observe the expected non-monotonic dependence of the 
exponent a on Pf Quantitatively there are deviations which arise from the fact that we 
only fit across three different values of c, and that (compare Fig. [T5l) we are not yet in 
the asymptotic low-c regime where the different bumps become well separated in time. 

We note briefly that in the paste-all regime, where the dynamics no longer separates 
into discrete stages, one can argue as above that the unperturbed defect density decays 
as n{t) ~ Perturbing T = ln(l/c) to T' = ln(l/c') = T + T^{(3h) + ... then 

shows that 

- XE{t, 0) ~ T2(lnt)t-^''^2 _ Tiyt2-''' oc Tutn{t). (124) 

At constant Ut this is proportional to T = ln(l/c): this logarithmic scaling is effectively 
a "smoothed" version across bumps of the dependence on with a = 1 — G [0, 1] in 
the irreversible coarsening regime. 

We can now combine xe and Ce into FD plots. Sample results from the numerical 
simulations are shown in Fig. [181 Three regimes can be discerned. Initially, for t^ 
close to t, XE(t,t^) is positive and of order c. This is a quasi-equilibrium response 
which we ignored in our estimates above. There follows a section where the predicted 
negative activated response becomes dominant; \xEit, tw)| increases here as t^ decreases, 
as expected from the dependence on 1 — t„/t. For small c, this factor can get close to 
unity before leaves the plateau that t is in. Only when this happens does CE(t,t^) 



Non- equilibrium dynamics of spin facilitated models 



48 



start to change significantly, and we enter the third regime where X-b(^, ^w) is essentially 
constant and the FD plot therefore horizontal. (Within this regime the FD plot then 
eventually reverses direction in the peculiar initial stage where is of order unity, when 
CE(t,t^) starts to increase again.) 

Our above theoretical estimates do not allow us to predict the precise behaviour of 
the FD plot around the initial quasi-equilibrium regime. Nevertheless, the fact that xe 
reaches values of order c suggests scaling both axes of the plot with c~^; Fig. [T8l indicates 
that a limit plot is then approached as c gets small at constant z/j. This describes the 
initial quasi-equilibrium regime and the crossover to negative values of xe- In the second 
regime that follows, the numerically obtained FD plots are straight with negative slopes 
of order unity. To rationalize this, consider t and t^ within the same plateau. Both 
are then small compared to the timescale governing the following stage of the dynamics 
{k = kt + 1), and we can linearize the evolution of the CE{t, tw) in both times. This then 
leads to a linear dependence on t — of ACE(t,t^) = CE(t,t) — CE(t,t^). Combined 
with the similarly linear relation xe oc (t — t^)/t, the FD plot should then be a straight 
line as observed numerically. Working out the relevant prefactors oit — t^ for AC e and 
Xe one can predict for the energy-FDR in this regime as (we omit the details) 

'1 



XE = -{kt + l) 




(125) 



Here the averages are over the domain size distribution, T{d) is the relaxation rate of 



domains of size d (see Appendix B ) and all quantities are evaluated at the beginning 
of the current plateau. The numerical value of (11251) can be found relatively easily for 
the first few plateaus. For kt = 0, there is only one nonzero rate, r(2). It follows 
that Td/r = 2, and evaluating the remaining moments d and rf^ of Po{d) gives the 
prediction Xe = —2.54. . ., in full agreement with the value Xe = —2.55(2) obtained 
from the simulations shown in Fig. [181 For kt = 1 one predicts similarly, by using that 
r(3) = 2c^/3 and r(4) = c^/4 (see [Appendix Bj ), the value Xe = —3.79 . . ., again in 



agreement with the simulation estimate Xe = —3.77(5). In the paste-all limit, where 
T{d) is significant only for the shortest domains, one has Td/T = d^i^ so that the energy 
FDR Xe = —{kt + 1)/[1 — exp(— 7)] grows linearly with the index of the stage of the 
dynamics. Intuitively this dependence arises because timescales grow as c~^^~^ and as 
kt grows so does the perturbation arising from the change in c. 

So far we have discussed the initial quasi-equilibrium regime of the energy FD plot 
and the straight line section with negative FDR Xe that follows. Because the "height" 
of the FD plot scales as X-b(^,0) ~ c^~"* and Xe is of order unity there, this section 
only extends by a small amount ~ c^~"* of the same order along the C^j-axis. In the 
remaining "third" section, Xe vanishes for c — as explained above. More precisely, 
this section can be divided further into subsections where the FDR scales as Xe ~ c, 
Xe ~ etc (up to c^'+^) as t^ moves backwards through the various plateaus of the 
dynamics. To see this, recall that Xe = —{dxE/dt^)/{dCE/dt^). The numerator 
equals -dxE/dt^ = XE{t,0)/t from f[T22[) which scales as c'^-'^'+C't+'^t) = ^i+kt^ ^^e 
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Figure 19. Non-local equal-time correlations Cr(t,t) from simulations (symbols) for 
T = 0.15 and different values of i/t; these compare well with the irreversible coarsening 
predictions shown as dashed lines. Inset: paste-all prediction for the same quantity, 
from Eq. ([^0]) . 



correlation function we do not have an explicit result that interpolates between plateaus 
but it seems reasonable to assume a smooth dependence on This then gives the 

estimate dCE/dt^ ~ {d / dt^)n{t^) which from flllSp will scale as c^""^^ to leading order. 
Putting this together gives Xe ~ c^^-^^ and hence scalings ~ c, c^, . . . , c^i+^ as claimed. 
The case where t„ is still in the same plateau as t {k^ = kt) is included and leads to the 
FDRs of order unity discussed above. 

3.3. Non-local observables 

In the previous subsections we found that both for local and global observables the 
correlation functions have values of order unity and exhibit a plateau structure that 
reflects the splitting of the dynamics into discrete stages (in the irreversible coarsening 
regime). The susceptibilities of local and global observables differ more strongly: for the 
local case, xo is always positive and of order c, and the local FDR is positive and of the 
order of c/n^{t^) (cf. the discussion at the end of Sec. 13.11) . The global susceptibility, 
on the other hand, is negative save for an initial quasi-equilibrium part of 0(c), and 
the associated FDR Xe is negative with values of order unity before dropping to zero 
- more precisely 0(c), O(c^), . . . - in the final segment of the FD plot. The question 
naturally arises of how these results relate to each other. 

Based on experience with the FA model we initially experimented with Fourier 
component observables. These, however, lead to FD plots that are very difficult to 
interpret; also, the limit of short wavevectors does not give direct access to the quantities 
for local observables. Instead we consider observables defined by random staggered fields 
ej with Gaussian correlations (ejej+r) = exp(— r^/2£^) in space. These have correlation 
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Figure 20. Non-local correlation function Crit^t^) /[n{t)n{t^)] for various negative 
r, plotted against v^. Symbols: simulation data for T = 0.15, for fixed I't = 1.2 (open 
symbols) and i't — 1-5 (full symbols). Dashed lines: theoretical prediction in terms of 
equal-time correlations, Cr{tw,tw)/n'^{tv!)i from eq. (|132p . 



and susceptibility 

Q(t,t„)= 5^e-^'/2^'a(t,tw), (126) 

r 

X.(t,tw) = $^e-'-'/2^'x.(t,tw). (127) 

r 

As I is increased from to oo we can thus directly interpolate between local and global 
observables. Before looking at this ^-dependence we turn to the non-local functions Cr 
and Xr as these are the building blocks from which and xt ^ire constructed. 

The equal-time correlations Cr{t,t) are simple to predict in the irreversible 
coarsening regime. Because the dynamics is of an independent interval type, domains 
from the domain size distribution Pk{d) for the current plateau k = kt are arranged in 
random order. Summing the probability of having an up-spin one, two (etc) domains 
along one gets 

Cr{t, t) = {mifym+rit)) - n\t) (128) 



n{t) 



-nit) + 5r,o + Pk{r) + V Pk{d)Pk{r - d) + . . 



d 



(129) 



This is simple to evaluate by numerical convolution given that we know Pk{d)] only the 
first few terms are needed since the m-th order convolution term is nonzero only for 
r > mdj^in{t) = m{2^ + 1). The result is shown in Fig. [19] for the first few plateaus 
and is in very good agreement with numerical data. The negative initial section for 
< r < dminit) arises because no up-spins exist in this distance range. This strong 
short-range repulsion also leads to oscillations in the correlation for larger values of 
r, including weak anti-correlations at some distances. In the paste-all limit one can 
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Figure 21. Scaling function H{x,9) for tlic non-local correlation function in the paste- 
all regime. Curves are for 6 — 1, 1.2, . . . , 3 from left to right. Lines: Simulations of 
paste-all dynamics; symbols: analytical prediction ()134p . 



exploit that the infinite sum of convolutions in f ll29p becomes a simple geometric series 
for the Laplace transforms. The result is proportional to exp(Ei(s)), and inverting the 
transform gives in terms of the functions (!98|) 

n{r/d^,^{t), 1), n{x, l) = -l + e^Yl j^M 



n^{t) 



i>i 



(130) 



A plot of this is included in Fig. [19] and clearly shows that small negative values occur 
e.g. around x = 2, outside the main "repulsion zone" < x < 1. 
The two-time correlations 

Cr(t,t^) = {ni+r(t)ni(t^)) -n{t)n{t^) (131) 

are more difficult but can still be obtained relatively simply for r < 0. The survival of an 
up-spin at site i + r from time t„ to time t depends only on the arrangement of domain 
lengths to the left of this site; in particular, it does not depend on whether an up-spin is 
present at site i (to the right of i-|-r) at time t^. Since the overall survival probability of 
an up-spin is just n{t)/n{t„), this implies (r;.j+r(^)^j(^w)) = [n{t)/n{t^)]{ni^ri'tw)ni{t^)) 
or, in terms of Cr(t,t^), 

a(Mw) ^ Cr{t t ) ^ ^ ^ 

n{t)n{t^) n\t^) ^ ' 

This prediction is very well verified by our simulation results in Fig. [20l which confirm 
for various negative r that the left hand side of Eq. fll32p is both independent of t and 
varies with in the matter predicted. 

For positive values of r, no similar simple argument applies. However, one can 
still calculate Cr{t,t^) within the paste-all regime; this involves keeping track of which 
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Figure 22. Non-local correlation function Cr{t,tv,), as obtained from simulations at 
ft — 1.5 and z^w — 0.5 and for the three temperatures indicated. The axes are as 
in Fig. [211 showing the correlations scaled by n{t)n{t^) against the scaled distance 
X — r / dnunitw) ■ The solid line graphs the paste- all prediction for the limit of large 
and Vf 



domain lengths from time have been swallowed up into the larger domains at time t 
Appendix E). One finds the following scaling: 



see 



X 



(133) 



n{t)n{t^) drninit^)' n(t) ' 

The scaling function Ti can be calculated explicitly in the range < x < 3 and 1 < < 3, 
with the result: 
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(134) 



(Notice that for negative distances, the scaling is independent of t from fll32p : 
TC{—x, 6) = H{—x, 1), and due to the spatial symmetry of the equal-time correlation this 
is identical to H{x, 1) from (11301) .) We plot H{x, 0) for positive x in Fig.[2T]for a range of 
values of 9, comparing also with data from direct simulations of the paste-all dynamics. 
The two-time correlations are seen to have an extremely rich spatial structure, which 
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Figure 23. Non-local susceptibility from simulations at T = 0.15 and for i/t = 0.7 
(top) and vt ~ 1.5 (bottom). We plot Xr{t,t^) against r for a range of waiting times 
as indicated by the legend. 



arises from the interplay of the different lengthscales dmin{t^) and d^i^(t). The strict 
exclusion zone for r < d^i^{t^) (x < 1) remains at the later time t. Larger distances r 
can separate spins at the two different times tw and t [Ti > — 1), but correlations only 
become positive beyond the lengthscale r = (ijnin(^) = Simulation data in the 
first plateaus (Fig. [22]) qualitatively follow these predicted trends; later plateaus would 
evidently need to be considered to see more quantitative agreement with the paste-all 
limit of large fc„ and kt. 

Next we consider the non-local susceptibility Xr{t,ty^). Here again the case of 
negative r is simpler: Xri^i tw) then has to vanish. This is because the evolution of spin 
Ui^rit) is only affected by the spins to its left; in particular, it remains unperturbed 
by a field applied at site i > i + r. For positive r, on the other hand, Xr(t,tw) will 
be nonzero but there are competing effects governing its sign and magnitude, making 
theoretical analysis difficult. We show selected simulation data for Xr{'t,t^) in Fig. [231 
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Figure 24. Scaling exponent obtained from fits of simulation data to Xr{t, 0) ~ c° 
at fixed Vf The exponent is plotted against i>t for r as indicated by the legend. 



and the exponent ar extracted by fits to the expected scahng Xr{t,0) ~ c"'' at fixed ut 
in Fig. El 

To develop some intuition, consider r = 1 and a local spin configuration 
{n^i,no,ni) = (1,0,1) at time tw = when the field is applied to site 0. This will 
increase the chances of uq fiipping up, after a time ~ c~^, and hence increase the 
probability of ni fiipping down. The resulting contribution to xi(^;0) is a "bump" 
around t = c~^ with negative sign and amplitude of order unity. If the next up-spin 
to the left of uq is further away, a similar effect occurs but with a longer delay because 
an up-spin front needs to propagate to Overall this should give rise to a series 

of negative bumps in Xi(^)O) centred around t = c~^, c^^ etc, similarly to XEit,0). 
At constant Ut the scaling with c is then ~ c^~°* as argued after (11231) . This is in 
reasonable agreement with the data (Fig. [241) . For r = 2, an argument based on initial 
configurations such as (n_i, no, ni, n2) = (1,0,0,1) suggests a similar structure but 
with the bump around t = c~^ absent, giving the scaling X2(^,0) ~ c^"*^* for z/t < 2 
(and ~ c^~"* thereafter). This is again in general accord with the numerical trends seen 
in Fig. EH Closer inspection reveals, however, that X2(i^,0) is in fact positive below 
Vt ~ 0.8. This sign change, which also makes the measurements of the exponent in 
this region rather unreliable, is not accounted for by our naive argument. For r = 3 
this issue becomes more pronounced still; in fact, Xal^jO) is positive throughout the 
range that we can explore (Fig. [23]) . The dominant contribution for t < c~^ is in this 
case an effect across domains: starting from (n_i, no, ni, n2, Tis) = (1,0,1,0,1), a field 
at site will speed up the disappearance of the domain formed by the first three spins. 
This causes ris to survive for longer^ giving a positive contribution to Xsl^^O). The 
timescale t ~ on which this term becomes significant is set by the rate for the first 
domain to disappear. One would expect, then, the scaling x^i^i 0) ~ c^~"* for ut < 1. In 
fact, the relevant scaling function - which for short chains can be calculated explicitly 




Aa(t,tw) 

Figure 25. Normalized FD plot for observables defined by Gaussian staggered fields, 
for T = 0.15 and vt = 0.5 (top), vt — 1.5 (bottom). The lengthscale I of the field 
correlation is given in the legend. 

- turns out to start off quadratically for small values of its argument, so that instead 
X3(t,0) ~ c^{i-at)^ in good agreement with our numerics (Fig. IMl) . It is not clear to 
us at present how to systematically account for effects of this type to get an overall 
prediction for the low-c scaling of Xrii, 0). Also the inclusion of the full tw-dependence 
of Xr(t,t^) is not trivial: for kt = (Fig. [23] top) a simple proportionality to (1 — t^/t) 
can be checked to describe the data well, whereas for kt = 1 (Fig. [23] bottom) this is no 
longer the case. 

Finally we show in Fig. [25] the FD plots for the ^-dependent correlation and 
response defined in f ll27p . These interpolate between the local {£ = 0) and global 
limits {£ oo) as they must. For the times considered here the dynamical lengthscales 
(typical domain sizes) are short enough that the global behaviour is approached already 
for moderate values of i. The approach to the local limit £ = 0, on the other hand, 
must be expected to become singular for small c. This can be seen as follows. Consider 
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Xe(^,^w) = X^rXr(^!^w)- We know that for = and fixed ut this quantity scales 
as c^~"* for c — > 0. There must therefore be one or several values of r for which 
Xr{t,0) scales in the same way. But then these values will dominate the sum over r 
in 0); in particular, they will always dominate over the local term (r = 0) with 
its 0(c) contribution. Thus, for any nonzero i we expect x^(^)O) ~ c^^"% whereas at 
£ — itself the scahng is ~ c. Assuming that similar arguments apply also to the 
two-time quantities x^(^) ^w); one expects the FD plots for sufficiently small c to switch 
effectively discontinuously from the local FD relation ior i — to behaviour dominated 
by non-local effects for £ > 0. 

4. Discussion and Conclusions 

We have studied two paradigmatic kinetically constrained models (KCMs) of glassy 
dynamics: the FA (Fredrickson-Andersen) and East models. Deploying a variety of 
analytical techniques and comparing with detailed numerical simulations, we have 
analysed in particular the correlation and response functions during the aging after 
a quench to low temperature, along with the resulting fiuctuation-dissipation ratios 
(FDRs). Local as well as global observables were considered, and Fourier mode and 
Gaussian staggered field observables respectively allowed us to interpolate between these 
two limiting cases. 

In the FA model, with its effective dynamics of diffusing and coagulating defects, 
a clear physical scenario emerges: the asymptotic FDR for well-separated times 
t ^ t^ is negative. Its precise value depends on whether one is above or below the 
critical dimension dc — 2, with X°° = — 3 for d > 2 and X°° — — 37r/(67r — 16) in 
d — 1. The underlying physics is, however, the same: the negative FDR arises from 
the activated nature of the aging dynamics. Where a temperature increase would, in 
equilibrium, increase the defect density n, here its main effect is to speed up the decrease 
of n{t) and so reduce its value. The asymptotic FDR can be determined from observables 
of any wavelength or characteristic Icngthscale, although for local observables this is very 
awkward because the interesting aging effects are buried underneath a dominant quasi- 
equilibrium signal. Much better suited is the global observable, i.e. the energy or total 
number of defects, which produces a fiuctuation-dissipation (FD) plot that is close (in 
d— 1) or exactly equal (in d > 2) to a straight fine of slope X°°. 

The East model has a more strongly cooperative behaviour than the FA model, 
and a correspondingly more subtle pattern of violations of the fiuctuation-dissipation 
theorem (FDT). For the local observable, the FDR Xo(t, t^) is always positive but small, 
of order c/n'^{tyj), where c is the (small) equilibrium defect concentration. In the paste- 
all limit of large domain sizes, which corresponds to a long-time limit taken within the 
aging regime, we find an intriguing similarity with mean-field predictions for spin glass 
dynamics: a continuous hierarchy of relaxation timescales leads to a curved FD plot 
that is effectively composed of a sequence of infinitesimal straight line segments. 

For global observables, on the other hand, also the East model displays the negative 
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FDRs characteristic of activated aging dynamics. For times t and t^ in the same 
"plateau" of the dynamics, the FDR has (negative) values of order unity which we 
can predict theoretically; as t^ becomes smaller it then drops to 0(c), O(c^) and so 
on. The contrast to the local observable can be traced to the different scaling with c, 
and different signs, of the distance dependent susceptibilities Xr- Observables probing 
intermediate lengthscales £ can interpolate between the local and global cases, although 
we argued that for c the non-local effects would dominate for any i > 0. 

The apparent decoupling of the local response from activation effects in the East 
model can be understood as follows. The response of a spin to a local field is governed 
by the history of its facilitating neighbour on the left. The field does not affect this 
"clock" because of the directed nature of the East model, and so one does not get any 
speed-up effects: the susceptibility only reflects the direct influence of the field and is 
positive. For spins at some distance r > from the site where the field is applied, one 
has the opposite scenario. These spins have a pure "speed-up" response: the spin at the 
field site will be up more often, thus accelerating the dynamics of all spins on the right. 

We already alluded in the introduction (Sec. II. 4p to the wider implications of our 
results. The non-trivial FDRs we have found in the aging regime reflect the growth of a 
purely dynamic lengthscale and cannot be related to the (trivial) equilibrium properties 
of KCMs. For most observables the FDRs are negative, precluding an interpretation in 
terms of an effective temperature. The negative sign nevertheless has a clear physical 
interpretation as arising from the activated nature of the aging dynamics. In the FA 
case, also the actual value of the asymptotic FDR (in the sense of widely separated 
times if: ^ t^) is robust among observables probing the entire range of lengthscales, 
from purely local behaviour to system-spanning global observables. In the East model 
the situation is more subtle, with the FD behaviour of local observables decoupled 
from activation effects. The latter do show up, however, in appropriate non-local and 
global observables. Because the East model has a hierarchy of relaxation timescales, the 
(negative) value of the FDR then also varies depending on which stage of the dynamics 
is being considered. 

To what extent will our results apply to other KCMs? One suspects that the 
strictly directed KCMs might form a somewhat separate group with regard to FD 
behaviour, because the effects of a local field can never propagate back to the site where 
it was applied. Qualitatively similar effects to those seen here for the East model would 
therefore be expected in, for example, the analogous three-dimensional model studied 
in Ref. [IQ]. The local FD plot shown in Fig. [26] demonstrates that this is indeed the 
case. In d = 3 there is again a hierarchy of timescales in the non-equilibrium dynamics, 
with the defect density relaxing in steps as in Fig. [61 This leads to FD plots that are 
(approximately) piecewise linear, with a number of segments depending on the position 
of t in the hierarchy (compare Fig. [10] for d = 1). The two-dimensional triangular 
plaquette model |126j also shows some broad similarities to the East model, but with 
additional subtleties because one now has to choose between studying the dynamics of 
the defects or of the underlying spin system. When the direction of facilitation is itself 
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Figure 26. FD plot for the local observable in the three-dimensional analogue of the 
East model, for T = 0.15 and i/f — 0.5, 1.5, 2.5 (from bottom to top). As in d = 1, 
the FD plot consists of approximately straight segments whose number increases with 
t (see Fig. [101). 

treated as a dynamical variable ^37j, the model becomes isotropic overall and one might 
expect the connection between local and global observables to be restored as for the FA 
model. It would therefore be very interesting to study the aging dynamics of such a 
model. 

A significant open question concerns the FD behaviour of undirected but highly 
cooperative KCMs, for example the FA model in d = 2 with two rather than one up- 
spin neighbours required to facilitate a move. Numerical studies on similar models with 
conserved dynamics |106l 1107] indicate a mean-field-like scenario without any apparent 
activation effects. It seems clear, however, that such effects should be present, certainly 
for global observables: the relaxation to lower defect densities must, as for our simpler 
one-spin facilitated FA models, proceed more quickly for higher c and result in a negative 
susceptibility. Here and also more widely in experimental studies we hope that our 
study will stimulate the search for well-defined negative FDRs. Arising as they do from 
activated dynamics, they should be essentially ubiquitous among systems exhibiting 
glassy dynamics. 
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Appendix A. Functions used in Section 12.21 

Two functions appearing throughout Sec. 12.21 are the modified Bessel function In(t) of 
integer order n \127\ 1146] , 

Ut)= r^^cos(ny,)e*^°^(^), (A.l) 

Jo ^TT 

and the function Hn(t), 

Hn{t) = I f dt'e-'' [In-lit') - In+lit')] . (A.2) 

Elementary properties that can be read off from these equations are /n(0) = 5nfl, 
I -nit) = Init), Hoit) = and H-nit) = —Hnit). In our analysis we make extensive use 
of the asymptotic scalings of these functions. At fixed order n and in the limit t ^ oo, 

e"*/„(t) ~ and Hn{t) ~ 1, (A.3) 

V 27rt 

whereas for t, n — ^ oo simultaneously with n^/t fixed, 

e-%(t)~-i=e-"^/(2*) and ~ (;^) • (A-4) 

Here $(2) = (2/y^) due^"^^ denotes the complementary error function as in the 
main text. The following functional relations are also of great value for our analysis, 

9i/„(t) = ^[/„_i(t) + 4+i(t)], (A.5) 

-I^[t)= -[In-lit) -In+lit)]. (A.6) 

To calculate Fourier transforms we use the identities (t > tw > 0) 

e-'^^In+min+jit) = e^5(^+i)'?/,_^. (2t cos I) , (A.7) 

n=— 00 
00 

^ ^ 6 Init t^)[In—m, ~l~ In+Tri\it ~\~ ^w) 
n=— 00 



A 

00 

e'"^''[In-l - In+l]it - tw)[/n-l " + t 



(A.9) 



where A = ^Jt'^ cos(g/2)2 + t^^ sin(g/2)2 and T„(x) = cos(n arccos x) is the Chebyshev 
polynomial of degree n. Equations flA.7p and flA.81) follow from the integral 
representation Eq. (lA.ip and trigonometric identities. The sum Eq. f lA.Qp can be reduced 
to Eq. (lA.SP by using Eq. (]A.6|) and expressing the resulting factor second 
derivative with respect to q. 
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Appendix B. Initial stages of irreversible coarsening dynamics 

In this section we give further details of the solution of the irreversible coarsening 
dynamics, leading in particular to the explicit expressions for the scaling functions flllQf - 
I12ip for stages A; = 0, 1, 2 of the dynamics. 

The dynamical equations ( IMl) for stage k can be solved in terms of the generating 
functions |116j 

G{z,t) = J2Pid,t)z', H{z,t)= Yl nd,t)z\ (B.l) 

d 2'''-i<d<2'= 

where H is the analogue of G restricted to the active domains. Bearing in mind that flMl) 
applies only for the inactive domain lengths d, one has 

^\G{z, t) - Hiz, t)] = Giz, t) (-^\ Hiz, t), (B.2) 



dt' ' ' ' ' ' ' ' ' \ dt^ 

and hence 

G{z, t) = l + [G{z, 0) - 1] exp[H{z, 0) - H{z, t)]. (B.3) 

In the limit t —* oo, where all active domains have disappeared {H — > 0), this gives a 
relation between the generating functions in plateaus k — 1 and k: 

Gk{z) = 1 + [Gk^iiz) - 1] exp[Hk^i{z)]. (B.4) 

If we define hk{d) as the inverse transform of exp[Hk{z)] — 1, i.e. exp[Hk{z)] — 1 = 
^dhk{d)z'^, this gives the following recursion for the domain length distributions 

Pkid) = Pk-M - hk^M + (Pfc_i * (B.5) 

Starting from the initial condition P^i{d) = 2^'^ one finds easily the particular values 
that we will need later: 

The overall distributions, which were obtained by evaluating the recursion numerically, 
are graphed in Fig. [51 

To get a prediction for the evolution of the defect density n{t) during stage /c, we 
exploit that l/n(t) = d{t) = {d/dz)G{z = l,t). From flRSl) this gives 

n{t) =n(0)exp[if(l,t) -/f(l,0)], (B.7) 

and we need the evolution of H{l,t). Since active domains cannot be recreated during 
the coarsening process, we can write 

H{l,t)= Yl Pid,t)= Y Pid,0)Sid,t), (B.8) 

where S{d,t) is the survival probability of a domain of length d. In stage k = this is 
very simple: the disappearance of a domain of length d = 1 corresponds to a mobile up- 
spin fiipping down, so S{l,t) = exp[— (1 — c)t] ~ exp(— t) for c ^ 0. Also -P-i(l) = 1/2 
at the beginning of stage 0, so H{l,t) = exp(— 1)/2. Inserting into (1B.7P directly leads 
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to the expression (11191) for (?o(C)- later stages and correspondingly larger d one can 
proceed using a "leaky" Markov chain that exits as soon as spin flips down. Taking 
= 1 as an example, active domains have length d = 2. The initial state of the domain 
is 101 (i.e. uq = 1, Hi = 0, = 1). The only other possible state that can be reached 
before spin n2 flips down is 111. The leaky transition matrix between these two states 
is 

•^'-(7 -2(r-e))- 

The flrst column represents the transition at rate c from 101 to 111; the second column 
has the reverse transition, with rate 1 — c, and in addition the rate 1 — c for exiting via 
111 — s> 110. The survival probability is just the probability of remaining in the chain, 
starting in state 101: 

S{2.t) = (1, l)e>« = + + (i - ^) e^", (B.IO) 




with D = a/1 — 2c + 5c2/4 and the eigenvalues Ai,2 = — 1 + c/2 ± D. To get gi{C) 
we need to consider the survival probability on the timescale t = (/c, taking the limit 
c — > at flxed (. Since A2 stays of order unity, the second exponential then disappears. 
In the flrst one, Ai = —c/2 to leading order while the prefactor tends to 1 so that overall 
S{2,t = (/c) — i> exp(— ^/2). Inserting into ( IB. 81) and using ( IB. 61) for the prefactor Po{2) 
then gives the result (I120p for gi{C). One proceeds similarly for (72(C) to derive (I12ip . 
the main difference being that there are now two active domain lengths, d = 3 and 
d = 4. The required leaky transition matrices are of size 4x4 and 8x8, respectively, 
and the survival probability is evaluated in the regime t = (/c^. By taking derivatives 
of the survival functions a.t ( = one further obtains the rates T{d) required for the 
prediction of the negative energy FDR (I125p . as 

r(i) = i, r(2) = |, r(3) = ^, m = j- (b.h) 

Appendix C. Exact relation between local correlation and response 

We outline the derivation of the exact relation ( 1 11 II) between the local correlation Cq 
and susceptibility Xo the East model. For the latter we have the expression (I106p . 
Xo(t,tvj) = c(l — c)[l — r(t,tvi)], in terms of the relaxation integral 

r(t,0 = (f(t,tj), f(t,t„) =exp (^-^ rft'n,_i(0^ . (C.l) 

For the correlation function we need to evaluate (ni(t)nj(t„)). As in the derivation of 
Xo we flrst take the history of spin nj_i from time to t as flxed. If also ni(t^) is given, 
the average value of Uiit) is c + [ni{t^) — c]f(t, t^), using the same argument as for xo- 
Thus 

(nj(t)ni(tw))n,{tw),n.-i(o...t) = ni{t^)[c+ (1 - c)f(t,t^)], (C.2) 
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where the subscript indicates the fixed quantities. Now we average over Uiit^), which has 
evolved from the initial average defect density n(0) at time to a later value governed 
by the relaxation integral f(t„, 0): 

(n,(t)n,(t^))„,_,(o...t) = {c+ [n(0) - c]f(t„, 0)}[c + (1 - c)f(t, t^)] (C.3) 
= + c(l - c)r(t, t^) + c[n(0) - c]f{t^, 0) 

+ (l-c)[n(0)-c]f(t,0), (C.4) 

where we have exploited that f{t,t^)f{t^,Q) = r(t,0). Still at fixed history of ni_i, the 
same arguments as above give the average densities at times t^ and t as 

(^i(tw))n,_i(0...i) = c + [n(0) - c]f{t^, 0), (C.5) 
(r^^(^))„,_,(o...^) =c+[r2(0)-c]f(t,0). (C.6) 
These can be used to eliminate the occurrences of f{t^, 0) and f{t, 0) from (10.4^ : 

{ni{t)ni{t^))n,_^^o...t) = + c(l - c)f(t, t^) + c[{ni{t„))n,_-,(o...t) - c] 

+ (l-c)[(n,(t))„^.,(o..,)-c]. (C.7) 

Averaging over the history of Ui-i gives then 

{ni{t)n,{t^)) = -cil - c)[l - r{t,t^)] + cn{t^) + {1 - c)n{t), (C.8) 
and so finally for the autocorrelation function Co(t,t^) = {ni{t)ni(t„)) — n{t)n{t^): 

Co(t, tw) = -c(l - c) [1 - r(t, tw)] + cn(tw) + (1 - c)n(t) - n(t)n(t„).(C.9) 

The first term on the r.h.s. is — xo(^;^w), and rearranging slightly gives the promised 
result (11 111) . It is clear from the derivation that this exact relation will hold for all 
kinetically constrained spin models with directed constraints, as long as we replace 
by the appropriate facilitation factor fi{t') for spin Ui. The key ingredient is 
that Ui does not affect the evolution of this facilitation factor; this is why we can first 
fix the facilitation history and average over the dynamics of n^. 



Appendix D. Energy correlations 

Here we sketch how the energy correlations for the East model can be calculated within 
the irreversible coarsening regime. The energy E = N/d is inversely proportional to the 
average domain length. Since the latter has fiuctuations 5d of over N'^^"^ we can write 

N{5d{t)6d{t^)) 

^-^'^'-^= d^it)d^it.) ■ (^-1) 

The fiuctuations 6d arise from the corresponding fiuctuations 6P{d) in the domain size 
distribution. At equal times, t^ = t, the latter are easy to obtain because of the 
independent interval nature of the dynamics: the actual arrangement of domain lengths 
at some given time can be obtained by repeatedly sampling domain lengths from the 
(average) domain length distribution P{d) and lining them up along the chain until the 
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total length is reached. By a relatively simple combinatorial calculation one then 
finds for large the following expression for the covariance of the 6P{d): 

N{6P{d)6P{d')) = d [P{d)6d,d' - P{d)P{d')] . (D.2) 

Summing over d and d', this implies {[J2d^Pi^)]'^) = as it must be because the 
distribution is always normalized. By multiplying with dd' first and then summing, on 
the other hand, one obtains the variance of Sd as d{d'^—d'^) and hence the expression (11161) 
for the equal-time energy correlations. 

For the two-time correlations one exploits that the recursion flB.5|) for the domain 
size distribution in the plateaus of the dynamics holds whatever the initial shape 
of this distribution. We can apply it, in particular, to a domain size distribution 
Pk^{d) + 6Pk„{d) in plateau perturbed by a small fluctuation. Linearization in the 
small quantities SPk„{d) ~ N^'^/'^ then tells us how this fluctuation propagates into the 
later plateau kt. We can write the outcome in the generic form 

5PkXd) = Y,Mk,kMd')SPk„{d'), (D.3) 

d' 

with an appropriate (inflnite) matrix Mj^^^^ depending on both the initial and final 
plateau. Inserting into (ID. II) gives then 

CE{t,t^) = dfjl^ dd"M^,,Md')N{5P,^{d')5P,M")) (D-4) 



d,d',d" 

= dk!d^' E dd"M,,,„ {d, d') [P,„ {d')5,.,,„ - P,^ {d')P,„ {d")] (D.5) 
d,d'4" 

= E dMk.^M, d')P,M){d' - 4J. (D.6) 

To evaluate the sum numerically without explicitly calculating and storing M^tk^, we 
perturb the domain size distribution in plateau k„ from Pk„{d) to Pk^{d) [1 -|- e{d — dk^)] 
with some small "field" e, and find the resulting small change (proportional to e) in the 
average domain length in plateau kt- Multiplying by the prefactor d^^dki then gives 
the results shown in Table [2l 

A closed form solution for the propagation of perturbations from t„ to t can be 
obtained in the paste-all regime of large domain lengths. We work with the normalized 
lengths X = d/dmmitw) as before; the clock variable is 6* = d^iriit) / dminitw) so that at 
"time" 6 there are no domains of length x < 9. The distribution Pe{x) of domain lengths 
at time 9 obeys the master equation 

^Pe{x) = -5{x - 9)Pg{9) + Pe{x - 9)Pe{9). (D.7) 

The first term on the right captures the disappearance of domains of length x = 9; these 
merge with their right neighbours into larger domains as represented by the second term. 
(Because x > 9 always, Pe{x) has a step discontinuity at x = 6'; the Pe{9) in the master 
equation is to be understood as the nonzero probability PqIx = 9^) to the right of this 
discontinuity, i.e. the probability density of the shortest domains present.) The initial 
condition Pi{x) at 6^ = 1 is the scaling distribution P{x) from (^71) . 
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To solve flD.7p one goes to Laplace transforms 

POO 

Peis)= / dxe-'^'Peix), (D.8) 
Jo 

to find 

^Pe{s)=e-'^Pe{e)[Pe{s)-l]. (D.9) 
Integrating w.r.t. 6 gives explicitly 

Peis) = 1 + [Pi{s) - I] exp d9'e-'''Pei9')^ . (D.IO) 

The initial condition is the Laplace transform of P{x), given by A(s) = 1 — exp[— Ei(s)]. 
Since the paste-all limit is a scaling regime one expects Pe{x) = Pi{x/9)/9 and in 
particular Pe{9) = Pi{l)/9 = 1/9. This is indeed the correct solution of (ID.lOP as the 
integral over 9' then becomes Ei(s) — Ei{9s) so that Pe{s) = 1 — exp[— Ei(^s)] = Pi{9s). 

To apply the above description of the paste-all dynamics to the calculation of the 
energy correlation function, we substitute Pg{x) Pe{x) + 5Po{x) everywhere and 
linearize in the small perturbation 5Pe. In the Laplace transform version ( ID. 101) this 
yields 

6Pe{s) = 5A(s)e^'(^)-^'('^) + [Pi{s) - iJe^iW-^He^) J d9' e-'''6Pe'{9'), (D.ll) 
or 

e'''^''^6Pe{s) = e^'^^'^SPiis) - J d9' e-'''6Pg>{9'). (D.12) 

As explained above, we need to insert as the initial perturbation 6Pi{x) = Pi{x){x—x) = 
Pi{x){x — e'^), i.e. 

SP^{s) = -£A(s) - e^A(s) = ^e-^'(^) - e^[l - e'^'^^)]. (D.13) 

This gives for the Laplace transform (1D.12P of the perturbed domain size distribution 
at time 9 

p-s re 

e^'(^^)5Pe(s) = - e^[e^'W - 1] - / d9' e~'''5Pe>{9'). (D.14) 

s Ji 

It remains to determine 6Pgr{9') from the condition that 6Pg{x) = for x < 9. On the 

l.h.s. of ( 1D.14P we have the (Laplace transform of) a convolution of Pg{x) with another 

function; this convolution then also vanishes for x < 9. The same must therefore be 

true of the r.h.s., which is the Laplace transform of 

e(x - 1) - ^ ^fi{x) - e{9 - x)5P,{x). (D.15) 

i>i 

The condition that this must vanish for x < 9 implies that in this range 

6P,ix) = Q{x - 1) - ^ ^ fi{x). (D.16) 

Z>1 
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Because we can make 6 as large as desired, this expression must then hold for all x. 
We can now reinsert this result into flD.14|) and take the limit s ^ to obtain the 
perturbation 6x0 in the average domain length, using that SPg{s) = —s6xg + 0(s^) and 
Ei(s) = -lns-7 + s + 0{s^): 

-^^ = -2 + e'- j\e'6Peie% (D.17) 
or, after inserting (1D.16P 

6xf) = e^6 



l^l-e^ + e-e^^jl' dd' fi{d')^ . (D.18) 



From (ID.6P we can then finally write down the two-time energy correlation. By 
considering the combination CE{t,t„)/n(t) we remove one factor of d^^ and also ensure 
that all factors of c/min(t„) arising from our length rescaling cancel from the result. The 
remainder of the prefactor is 1/ {e"'e"'9) so that overall 

CE{t,t^)/n{t) = e-\e + l)-l-Y,\ f de'fiiO'). (D.19) 

This is the scaling function Q{9) given in ( lllTP ; for < 2 it can be written explicitly as 
g{e) = e-T(0 + 1) - 1 - In e. (D.20) 

For ^ = 1 in particular one obtains the scaling of the equal-time correlations in the paste- 
all limit as CE{t, t)/n{t) = 2e~'^ — 1 = 0.1229 . . . This value gives the y-axis intercept of 
the paste-all curve in Fig. [Ti 



Appendix E. Non-local correlations 

In this appendix we describe how the non-local two-time correlations in the East model 
can be calculated in the paste-all regime. The aim is to derive the scaling function 
T-C{x, 6) defined in (I133p . in the nontrivial region of positive x. 

As explained in the main text, two-time spatial correlations require one to keep 
track of which domains from time t^ have been merged into the domains at time t. 
We can then characterize a domain at t by the (scaled) lengths of the t^-domains it 
contains. If these lengths are, in order, xi, x;, we write the fraction of such t- 
domains as Pq\xi, . . . ,x/); here 6 = n{t„)/n{t) = (imm(i)/'^mm(^w) as before. For this 
more detailed description of the domains, the master equation ( ID. 71) becomes 

^Pf (xi, ...,xi)= - 6{x^ + ... + XI- ^)Pf (xi, ...,xi) (E.l) 

i-i 

+ '^6{xi + ... + xi' - e)P'^^ \xi, . . . ,Xi/)Pj'"' \xv+i, . . . ,xz). 
i'=i 

The first term represents the usual disappearance of domains of total length 9. The 
second one keeps track of the appearance of new domains; as two domains merge, the 
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lists of the lengths of the t„-domains which they contain are simply appended to each 
other. 

It turns out that the above equations for the Pq^ can be solved explicitly, but 
it is useful to check first how they enter the scaling function T-C{x,6). From the 
definition (11331) . H + 1 is the density of tw-spins a distance —x from a t-spin, divided 
by the overall density of t^-spins. The latter is 1/x = in the paste-all regime, so 
e~'^(7i + 1) for positive x is the density of t^-spins that are a distance x to the left of a 
i-spin. This can be written as 

e-\n + 1) = Qe{x) + {Qe * Pe){x) + (Q, * * P,)(x) + . . . (E.2) 

The terms in this series represent tw-spins that are within (or on the left boundary of) 
the first, second, third etc. t-domain to the left of the t-spin considered, hence the 
appearance of the length distribution PqIx) of t-domains. The factor Qe{x) records the 
positions of t„-spins within t-domains: 

oo 

Qe{x) = Y.Qt^\x), (E.3) 

m=l 
oo „ 

Qe"\x)= / dxi---dxiPg\xi,...,xi)5{xi-rn+i + ■■■ + xi-x),{EA) 

l=m 

where Q^^^ (x) accounts for the contributions from the m-th t„-spin within a t-domain, 
counting leftwards from the t-spin on its right boundary. We see that we do not need 
all details of the Pg'''; the Qg"*^ are sufficient. Unfortunately, the latter do not obey 
closed equations because they do not contain information about the total length of 
the t-domain. We therefore generalize so that this is also kept track of and define the 
quantities 



oo „ 

Q^^\x,Xtot) = ^ i dxi---dxiPf\xi,...,xi) 

l=m 



X 



X (5(x;_m+l + ■ ■ ■ + Xi - X)5{xi + . . . + Xi - Xtot)- (E.5) 

These give the fraction of t-domains of length Xtot for which the m rightmost t„-domains 
contained within add up to a length x. Taking the simplest case m = 1 as an example, 
one now derives from ( lE.ll) the evolution equation 

d ~ ~ 1 ~ 

■^Q^g\x,Xtot) = - S{xtot - 6)Q^\x,Xtot) + -Qf\x,Xtot " 6)- (E.6) 

The factor 1/9 here is the total rate of disappearance of domains; formally it is calculated 
from 



oo „ 

/ dxi ■ ■ ■ dxi' 6{xi + . . . + xi> — 6)Pg \ 
i'=i ^ 



xi,...,xi'), (E.7) 



which just gives the density Pe{d) = 1/0 of t-domains of length 6 (as enforced by the 
delta function) at time 6. 
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The initial condition for ( ]E.6l) at 6* = 1, i.e. t = tw, is Q^i\x,Xtot) = Pi{x)S{x — 
Xtot) = P{x)6{x — Xtot) because there is then no distinction between tw-domains and t- 
domains. (It is only through this initial condition that the x-dependence enters.) Going 
to Laplace transforms Xtot s, integrating with respect to 6 and transforming back 
gives then 



e ^ 



Xtot ~ X 

e 



Qg '* {x, X ) — P{x)S{XiQt — x) 



(E.8) 



Xtot 



9' 



Qi\x,0') 



where 



£{x) = Six) + J2lfii^) 
i>i 



(E.9) 



is the inverse Laplace transform of exp[Ei(s)]. We now need to find Q^fl,\x, 9') from the 
condition that Q'f\x^x') vanishes for x' < 9. Since £{■) is zero for negative argument, 
the l.h.s. of (lE.SP vanishes for xtot < so that in this regime 

9'^ 



P{x)£{xt. 



,ot 



XI 



d9' 



Xtot 



9' 



q£^(x,^')- 



(E.IO) 



Again because of the vanishing of £{■) for negative argument we can restrict the 
integration to 9' < Xtot- Also Qlf}{x,9') must vanish for 9' < x - the total length 
cannot be smaller than the length of the rightmost tw-domain within - so the lower 
integration limit can be raised from ^' = 1 to ^' = x. After relabelling Xtot 9 our 
condition becomes 



P{x)£{9~x) 



d9'f9-9' 



9' 



-£ 



9' 



g2^(x,0') 



(E.ll) 



It is from this expression that Q^q) {x^ 9') is to be found. The simplest regime is 9 /x < 2. 
Then the argument of £{■) on the right is always less than 1 and so only the first term 
in ( IE. 91) contributes, giving 

gW(a;, 9) = P{x)£{9 - x), for 9/x < 2. (E.12) 

Since x > 1 this solution applies, in particular, whenever 9 < 2, where it simplifies 
further to 

Q^g\x, 9) = P{x)5{9 - x). (E.13) 

This makes sense: a t-domain of size 9 < 2 cannot contain two or more tw-domains, so 
that necessarily 9 = x. Inserting (lE.lSP into (IE. 80 yields 



^dx^ 
9 ^ 



Xtot 



X 



9 



Q^^\x,x') = P{x) 



£{Xtot — x) — 



Q{9 



x] 



-£ 



X 



Xtot 



X 



X 



(E.14) 



One can now integrate over Xtot, or equivalently go back to Laplace transforms and take 
s — s> 0, to find 



dx'Q^^\x,x') = 9P{x) 



e( 



X] 



X 



(E.15) 
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This expression is exact for 6 < 2. If also x < 2, then it is the only contribution to (IE. 20 
because all other terms involve two or more t„-domains whose combined length will be 
above 2. We have thus obtained the first line of f ll34p . To get the full expression valid 
for ^ < 3 and x < 3 one has to extend the above calculation for Q^f^\x) to the range 
2 < 9 < 3] one also needs to calculate Qg (x). We omit the details. 



Appendix F. Persistence function 

Here we outline how to calculate the persistence function of down-spins f lllOp in the 
paste-all regime of the East model. From the discussion before fllOQp this governs the 
behaviour of the local susceptibility xo the East model. 

As discussed in Sec. |3l in the paste-all regime each domain contains an equilibration 
zone on its left, of size 9. (We use the same definitions of x and 9 as in previous 
sections.) Each equilibration zone has been fully "swept" by up-spins, i.e. it contains no 
persistent down-spins. However, regions to the right of the current equilibration zone 
may already have been swept previously. For example, when two domains merge, only 
the equilibration zone of the original domain on the left remains, but we have to keep 
track of the fact that the equilibration zone in the domain on the right has already been 
swept by up-spins. We therefore define Pg{u,x) as the fraction of domains that have 
length X and an unswept zone of length u (on their right end). This has two components: 

Pe{u, x) = Pe{x)6{x - 9 - u) + P^'(m, x). (F.l) 

The first part represents domains where the swept area is exactly the equilibration zone 
and so u = X ~ 9; the second part describes domains where the swept zone is larger 
and correspondingly u < x — 9. Initially, at ^ = 1, only the first part is present, and 
so P{{x) = P{x) and P['{u,x) = 0. This is because we want to count persistence from 
time tw, i-e. the swept zones are reset to the current equilibration zones at that point. 
At some later "time" 9, the fraction of the chain occupied by persistent down-spins is 
{u)0 divided by the average domain length e'^9. The function !F{9) is the complement 
of this, i.e. the fraction of swept areas 

m = 1 - (F-2) 

To calculate the required average {u)o we start from the evolution equations for P' 
and P", 

^^P',{x) = -6{x- 9)P',{x) + P;(x - ^, x), (F.3) 
^p;(m,x)= -6{x-9-u)P'^{x-9,x) + ^-Pe{u,x-9). (F.4) 



In (IF.3P the first term describes disappearance of domains of length 9\ note that 
only domains in the P'-part of the distribution can disappear in this way because at 
disappearance the equilibration zone covers the whole of the domain and is therefore 
identical to the swept zone. The second term accounts for the fact that when the 
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equilibration zone (length 6) grows long enough to cover all of the existing swept zone 
(length X — u = X — {x — 9) = 6), domains are counted in P' rather than P". The first 
term in (1F.4I) records the corresponding loss of domains from P". The last term, finally, 
accounts for the creation of new domains as domains of length 6 disappear and merge 
with their right neighbours. The factor 1/^ is the total rate of disappearance of domains 
as in flE.6l) : the length u of the unswept zone remains unchanged during a merger while 
the overall domain lengths 6 and x — 6 add. 

The first evolution equation ( IF.SP can be integrated directly to give 



Pl}{x) = Q{x - 



(F.5) 



p{x) + J do' p'^,{x-e\x) 

while for the second one we insert (IF. II) and then proceed as in the derivation of (lE.Sp 
to obtain 



dx' ^ f X — X 



Pf,{u,x) = I — 



X 



u 



Pa'iu.e' 



U 



1 ^fx-29' -u\ ^, 
+ 0. ]Pe'iO' + u) 



■ (F.6) 



One sees that all quantities are determined once we know Pg,{u, 6' + u). Fortunately, this 
vanishes for 6' < 2. The reason is that domains that have not yet merged have identical 
swept and equilibration zones and so are counted in P'. Domains in the P"-part of 
the distribution must then be the product of a merger of at least two of the domains 
that were present at time Their swept zone is, as a consequence, the result of a 
merger of at least two equilibration zones and must have length at least 2. It follows 
that Pq,{u,x) = for u > X — 2, and substituting x = 9' + u with 6' < 2 proves our 
claim. 

With the above simplification for < 2 we get Pq{x) = 0(x — 6)P{x) and 



9 



or, after integration over x, 



d9' 



dO' fx-2e' -u 



dx'p'^{u,x') = e ^^p{e' + u). 

Putting everything together gives for the average of u 



^Pie' + u),iF.7) 



(F.8) 



dx duu Pg{u^x) 



dx ix - e)P'Jx) 



(F.9) 



dxduuPg{u,x) (F.IO) 
de' ^ 

1 



dx{x-9)Q{x-e)P{x)+ duu9 -—P{9' + u) (F.ll) 

J Ji (^') 

(^^(^'), (F.12) 
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where 

I{9) = I dx{x - e)P{x) (F.13) 



/oo 
dx{x - e)P{x) - j dx{x - e)P{x) (F.14) 

-9- [ dx (l--^ (F.15) 



J 1 \ X / 

= e'^ + i-2e + e\ne. (p.ie) 

Carrying out the remaining ^'-integration in (1F.12P yields 

{u)e = (e^ -l)9-9\n9 + ^ln'^e, (F.17) 

and, after inserting into (lF.2p . the result (11 101) from the main text. One can push the 
calculation further, certainly up to < 3, but the resulting expressions become very 
unwieldy and so are not given here. 
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